From 8f616867c83ffd1edabd392ae90940ddc86f1482 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Mon, 8 Jul 2024 09:14:53 -0400 Subject: [PATCH 01/25] Start adding NORDIC options. --- fmriprep/cli/parser.py | 7 +++++++ fmriprep/config.py | 2 ++ 2 files changed, 9 insertions(+) diff --git a/fmriprep/cli/parser.py b/fmriprep/cli/parser.py index 63b4cd9d..4c5c0aa1 100644 --- a/fmriprep/cli/parser.py +++ b/fmriprep/cli/parser.py @@ -447,6 +447,13 @@ def _slice_time_ref(value, parser): 'It is faster and less memory intensive, but may be less accurate.' ), ) + g_conf.add_argument( + '--use-nordic', + action='store_true', + dest='use_nordic', + default=None, + help='Apply NORDIC denoising to the BOLD data', + ) 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 6ba0ae1c..8365772e 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""" + use_nordic = None + """Apply NORDIC denoising to the BOLD data.""" class loggers: From 5bae8340f36a3e9fa67ed801bdfd38b9082fd3a9 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 7 Aug 2024 17:48:51 -0400 Subject: [PATCH 02/25] Keep working on NORDIC. --- fmriprep/workflows/bold/fit.py | 39 +++++++- fmriprep/workflows/bold/nordic.py | 159 ++++++++++++++++++++++++++++++ 2 files changed, 196 insertions(+), 2 deletions(-) create mode 100644 fmriprep/workflows/bold/nordic.py diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index f66d49eb..63022964 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -46,6 +46,7 @@ # BOLD workflows from .hmc import init_bold_hmc_wf +from .nordic import init_bold_nordic_wf from .outputs import ( init_ds_boldmask_wf, init_ds_boldref_wf, @@ -755,6 +756,27 @@ def init_bold_native_wf( 'Multi-echo processing requires at least three different echos (found two).' ) + if config.workflow.use_nordic: + # Look for (1) phase data, (2) magnitude noRF data, (3) phase noRF data + bids_filters = config.execution.get().get('bids_filters', {}) + + norf_files = get_norfs( + bold_series, + 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) + run_stc = bool(metadata.get('SliceTiming')) and 'slicetiming' not in config.workflow.ignore workflow = pe.Workflow(name=name) @@ -794,6 +816,10 @@ def init_bold_native_wf( ) outputnode.inputs.metadata = metadata + nordicbuffer = pe.Node( + niu.IdentityInterface(fields=['bold_file']), + name='nordicbuffer', + ) boldbuffer = pe.Node( niu.IdentityInterface(fields=['bold_file', 'ro_time', 'pe_dir']), name='boldbuffer' ) @@ -814,16 +840,25 @@ def init_bold_native_wf( (bold_source, validate_bold, [('out', 'in_file')]), ]) # fmt:skip + if config.workflow.use_nordic: + nordic_wf = init_bold_nordic_wf(norf_files, mem_gb=mem_gb) + workflow.connect([ + (validate_bold, nordic_wf, [('out_file', 'inputnode.bold_file')]), + (nordic_wf, nordicbuffer, [('outputnode.bold_file', 'bold_file')]), + ]) # fmt:skip + else: + workflow.connect([(validate_bold, nordicbuffer, [('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')]), + (nordicbuffer, 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')])]) + workflow.connect([(nordicbuffer, 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 new file mode 100644 index 00000000..002f85bc --- /dev/null +++ b/fmriprep/workflows/bold/nordic.py @@ -0,0 +1,159 @@ +# 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 + +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 + ------- + nordic_file + NORDIC-denoised BOLD 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=['bold_file', 'norf_file', 'phase_file', 'phase_norf_file'], + ), + name='inputnode', + ) + outputnode = pe.Node(niu.IdentityInterface(fields=['nordic_file']), name='outputnode') + + # Add noRF file to end of bold_file if available + add_noise = pe.Node(niu.IdentityInterface(fields=['bold_file']), name='add_noise') + + if phase: + # Do the same for the phase data if available + add_phase_noise = pe.Node( + niu.IdentityInterface(fields=['bold_file', 'norf_file']), + name='add_phase_noise', + ) + validate_complex = pe.Node( + niu.IdentityInterface( + fields=['mag_file', 'phase_file', 'n_mag_noise_volumes', 'n_phase_noise_volumes'], + ), + name='validate_complex', + ) + + # 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( + niu.IdentityInterface(fields=['bold_file', 'n_noise_volumes']), + name='split_noise', + ) + + workflow.connect([ + (inputnode, add_noise, [('bold_file', 'in_file')]), + (nordic, split_noise, [('mag_file', 'bold_file')]), + (split_noise, outputnode, [('bold_file', 'nordic_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')]), + ]) # 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')]), + ]) + + return workflow From df8ae03aeebef051e646cda80a51905353f0ca65 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 7 Aug 2024 19:38:19 -0400 Subject: [PATCH 03/25] Looks good minus actual NORDIC node. --- fmriprep/interfaces/nordic.py | 166 ++++++++++++++++++++++++++++++ fmriprep/workflows/bold/nordic.py | 38 +++---- 2 files changed, 185 insertions(+), 19 deletions(-) create mode 100644 fmriprep/interfaces/nordic.py diff --git a/fmriprep/interfaces/nordic.py b/fmriprep/interfaces/nordic.py new file mode 100644 index 00000000..5510f451 --- /dev/null +++ b/fmriprep/interfaces/nordic.py @@ -0,0 +1,166 @@ +"""BIDS-related interfaces.""" + +from pathlib import Path + +from nipype.interfaces.base import ( + BaseInterfaceInputSpec, + File, + SimpleInterface, + TraitedSpec, + isdefined, + traits, +) + + +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 diff --git a/fmriprep/workflows/bold/nordic.py b/fmriprep/workflows/bold/nordic.py index 002f85bc..9a03cf9f 100644 --- a/fmriprep/workflows/bold/nordic.py +++ b/fmriprep/workflows/bold/nordic.py @@ -32,6 +32,7 @@ from nipype.pipeline import engine as pe from ... import config +from ...interfaces.nordic import AppendNoise, RemoveNoise, ValidateComplex LOGGER = config.loggers.workflow @@ -74,8 +75,10 @@ def init_bold_nordic_wf( Outputs ------- - nordic_file + mag_file NORDIC-denoised BOLD series NIfTI file + phase_file + NORDIC-denoised phase series NIfTI file """ from niworkflows.engine.workflows import LiterateWorkflow as Workflow @@ -85,27 +88,25 @@ def init_bold_nordic_wf( """ inputnode = pe.Node( niu.IdentityInterface( - fields=['bold_file', 'norf_file', 'phase_file', 'phase_norf_file'], + fields=['mag_file', 'norf_file', 'phase_file', 'phase_norf_file'], ), name='inputnode', ) - outputnode = pe.Node(niu.IdentityInterface(fields=['nordic_file']), name='outputnode') + 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(niu.IdentityInterface(fields=['bold_file']), name='add_noise') + add_noise = pe.Node(AppendNoise(), name='add_noise') if phase: # Do the same for the phase data if available - add_phase_noise = pe.Node( - niu.IdentityInterface(fields=['bold_file', 'norf_file']), - name='add_phase_noise', - ) - validate_complex = pe.Node( - niu.IdentityInterface( - fields=['mag_file', 'phase_file', 'n_mag_noise_volumes', 'n_phase_noise_volumes'], - ), - name='validate_complex', - ) + 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( @@ -115,10 +116,7 @@ def init_bold_nordic_wf( ) # Split noise volumes out of denoised bold_file - split_noise = pe.Node( - niu.IdentityInterface(fields=['bold_file', 'n_noise_volumes']), - name='split_noise', - ) + split_noise = pe.Node(RemoveNoise(), name='split_noise') workflow.connect([ (inputnode, add_noise, [('bold_file', 'in_file')]), @@ -146,6 +144,8 @@ def init_bold_nordic_wf( ('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([ @@ -154,6 +154,6 @@ def init_bold_nordic_wf( ('n_noise_volumes', 'n_noise_volumes'), ]), (add_noise, split_noise, [('n_noise_volumes', 'n_noise_volumes')]), - ]) + ]) # fmt:skip return workflow From 1c7bed51984edb962a9987c6b3dcd83356826899 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 7 Aug 2024 19:52:55 -0400 Subject: [PATCH 04/25] More work. --- fmriprep/workflows/bold/fit.py | 70 +++++++++++++++++++++++++++++-- fmriprep/workflows/bold/nordic.py | 7 +++- 2 files changed, 71 insertions(+), 6 deletions(-) diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index 63022964..d06b56b2 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -760,8 +760,9 @@ def init_bold_native_wf( # Look for (1) phase data, (2) magnitude noRF data, (3) phase noRF data bids_filters = config.execution.get().get('bids_filters', {}) - norf_files = get_norfs( + norf_files = get_associated( bold_series, + suffix='noRF', entity_overrides=bids_filters.get('norf', {}), layout=layout, ) @@ -777,6 +778,43 @@ def init_bold_native_wf( ) config.loggers.workflow.info(norf_msg) + phase_files = get_associated( + bold_series, + suffix='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 = get_associated( + bold_series, + 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) @@ -841,10 +879,34 @@ def init_bold_native_wf( ]) # fmt:skip if config.workflow.use_nordic: - nordic_wf = init_bold_nordic_wf(norf_files, mem_gb=mem_gb) + 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')]), + ]) # fmt:skip + + 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')]), + ]) # fmt:skip + + 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')]), + ]) # fmt:skip + + nordic_wf = init_bold_nordic_wf(phase=has_phase, mem_gb=mem_gb) workflow.connect([ - (validate_bold, nordic_wf, [('out_file', 'inputnode.bold_file')]), - (nordic_wf, nordicbuffer, [('outputnode.bold_file', 'bold_file')]), + (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')]), ]) # fmt:skip else: workflow.connect([(validate_bold, nordicbuffer, [('out_file', 'bold_file')])]) diff --git a/fmriprep/workflows/bold/nordic.py b/fmriprep/workflows/bold/nordic.py index 9a03cf9f..c820db28 100644 --- a/fmriprep/workflows/bold/nordic.py +++ b/fmriprep/workflows/bold/nordic.py @@ -119,9 +119,12 @@ def init_bold_nordic_wf( split_noise = pe.Node(RemoveNoise(), name='split_noise') workflow.connect([ - (inputnode, add_noise, [('bold_file', 'in_file')]), + (inputnode, add_noise, [ + ('mag_file', 'bold_file'), + ('norf_file', 'norf_file'), + ]), (nordic, split_noise, [('mag_file', 'bold_file')]), - (split_noise, outputnode, [('bold_file', 'nordic_file')]), + (split_noise, outputnode, [('bold_file', 'mag_file')]), ]) # fmt:skip if phase: From fc3729ae94a90dcad1b16f188dcc3685b089a793 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 7 Aug 2024 21:05:51 -0400 Subject: [PATCH 05/25] Add get_associated function. --- fmriprep/utils/bids.py | 15 +++++++++++++++ fmriprep/workflows/bold/fit.py | 8 ++++---- 2 files changed, 19 insertions(+), 4 deletions(-) diff --git a/fmriprep/utils/bids.py b/fmriprep/utils/bids.py index 0bdc03ff..9363196c 100644 --- a/fmriprep/utils/bids.py +++ b/fmriprep/utils/bids.py @@ -414,3 +414,18 @@ 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.append(layout.get_nearest(source_file, strict=True,**query)) + + 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/fit.py b/fmriprep/workflows/bold/fit.py index d06b56b2..ca12c584 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -41,7 +41,7 @@ ReconstructFieldmap, ResampleSeries, ) -from ...utils.bids import extract_entities +from ...utils.bids import extract_entities, get_associated from ...utils.misc import estimate_bold_mem_usage # BOLD workflows @@ -762,7 +762,7 @@ def init_bold_native_wf( norf_files = get_associated( bold_series, - suffix='noRF', + query={'suffix': 'noRF'}, entity_overrides=bids_filters.get('norf', {}), layout=layout, ) @@ -780,7 +780,7 @@ def init_bold_native_wf( phase_files = get_associated( bold_series, - suffix='phase', + query={'part': 'phase'}, entity_overrides=bids_filters.get('phase', {}), layout=layout, ) @@ -799,7 +799,7 @@ def init_bold_native_wf( phase_norf_files = get_associated( bold_series, - suffix='noRF', + query={'part': 'phase', 'suffix': 'noRF'}, entity_overrides=bids_filters.get('norf', {}), layout=layout, ) From a21d415d88232dbbd4a0ed3ada4373be45f87b25 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 8 Aug 2024 15:39:14 -0400 Subject: [PATCH 06/25] Update nordic.py --- fmriprep/interfaces/nordic.py | 94 ++++++++++++++++++++++++++++++++++- 1 file changed, 93 insertions(+), 1 deletion(-) diff --git a/fmriprep/interfaces/nordic.py b/fmriprep/interfaces/nordic.py index 5510f451..de8043b9 100644 --- a/fmriprep/interfaces/nordic.py +++ b/fmriprep/interfaces/nordic.py @@ -1,8 +1,11 @@ -"""BIDS-related interfaces.""" +"""NORDIC-related interfaces.""" +import os from pathlib import Path +from string import Template from nipype.interfaces.base import ( + BaseInterface, BaseInterfaceInputSpec, File, SimpleInterface, @@ -10,6 +13,7 @@ isdefined, traits, ) +from nipype.interfaces.matlab import MatlabCommand class _ValidateComplexInputSpec(BaseInterfaceInputSpec): @@ -164,3 +168,91 @@ def _run_interface(self, runtime): 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 From ae127093e4d3962354cc2f507cb17a3444b427cc Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sat, 9 Nov 2024 09:14:35 -0500 Subject: [PATCH 07/25] 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 From 39c6c85db8603336563d29c123c2d81fecb9317d Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sat, 9 Nov 2024 11:19:23 -0500 Subject: [PATCH 08/25] Add mrtrix3 to environment. --- env.yml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/env.yml b/env.yml index d680f3a6..0be66d83 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 From c92f653b14cf2e4bea77026ef874741ab0d970de Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sat, 9 Nov 2024 11:28:38 -0500 Subject: [PATCH 09/25] Run ruff. --- fmriprep/utils/bids.py | 2 +- fmriprep/workflows/bold/fit.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/fmriprep/utils/bids.py b/fmriprep/utils/bids.py index 06219f6d..24d886e6 100644 --- a/fmriprep/utils/bids.py +++ b/fmriprep/utils/bids.py @@ -428,7 +428,7 @@ def get_associated(source_files, query, entity_overrides, layout): query.update(entity_overrides) associated = [] for source_file in source_files: - associated.append(layout.get_nearest(source_file, strict=True,**query)) + associated.append(layout.get_nearest(source_file, strict=True, **query)) if len(associated) not in (0, len(source_files)): raise ValueError( diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index f2804b90..01bff51d 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -43,10 +43,10 @@ ) 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 -from .denoise import init_bold_dwidenoise_wf from .outputs import ( init_ds_boldmask_wf, init_ds_boldref_wf, From 0e88b6cf6da2d4a7e93e3976112d1fedf7d4c6dc Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sat, 9 Nov 2024 11:35:06 -0500 Subject: [PATCH 10/25] Fix connection. --- fmriprep/workflows/bold/fit.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index 01bff51d..a4904f40 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -940,11 +940,11 @@ def init_bold_native_wf( bold_stc_wf = init_bold_stc_wf(metadata=metadata, mem_gb=mem_gb) workflow.connect([ (inputnode, bold_stc_wf, [('dummy_scans', 'inputnode.skip_vols')]), - (denoisebuffer, 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([(denoisebuffer, boldbuffer, [('out_file', 'bold_file')])]) + workflow.connect([(denoisebuffer, boldbuffer, [('bold_file', 'bold_file')])]) # Prepare fieldmap metadata if fieldmap_id: From 0f53b42fa4ac5b90630a22e424139e3d8082edeb Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sat, 9 Nov 2024 12:02:14 -0500 Subject: [PATCH 11/25] Test the new behavior. --- .circleci/config.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.circleci/config.yml b/.circleci/config.yml index bc2689be..19fa7819 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -712,6 +712,7 @@ jobs: /tmp/data/${DATASET} /tmp/${DATASET}/fmriprep participant \ --fs-no-reconall --sloppy --write-graph \ --output-spaces MNI152NLin2009cAsym \ + --thermal-denoise-method mppca \ --mem-mb 14336 --nthreads 4 --anat-only -vv --notrack fi - run: From 3c1f61a55cbce24f78514c65c1aa32c3dc621959 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sat, 9 Nov 2024 13:11:33 -0500 Subject: [PATCH 12/25] I passed it to the anat-only run. :( --- .circleci/config.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.circleci/config.yml b/.circleci/config.yml index 19fa7819..456fc375 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -712,7 +712,6 @@ jobs: /tmp/data/${DATASET} /tmp/${DATASET}/fmriprep participant \ --fs-no-reconall --sloppy --write-graph \ --output-spaces MNI152NLin2009cAsym \ - --thermal-denoise-method mppca \ --mem-mb 14336 --nthreads 4 --anat-only -vv --notrack fi - run: @@ -756,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: From c59d589fd53d42687c374f15aa16c0866278e173 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sat, 9 Nov 2024 13:31:40 -0500 Subject: [PATCH 13/25] Fix bug. --- fmriprep/utils/bids.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/fmriprep/utils/bids.py b/fmriprep/utils/bids.py index 24d886e6..2d98dd7e 100644 --- a/fmriprep/utils/bids.py +++ b/fmriprep/utils/bids.py @@ -428,7 +428,9 @@ def get_associated(source_files, query, entity_overrides, layout): query.update(entity_overrides) associated = [] for source_file in source_files: - associated.append(layout.get_nearest(source_file, strict=True, **query)) + 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( From c9f655d0c9b5994889d08ff4d64cf482f4849f47 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sat, 9 Nov 2024 16:03:40 -0500 Subject: [PATCH 14/25] Fix order of connections. --- fmriprep/workflows/bold/fit.py | 54 +++++++++++++++++++--------------- 1 file changed, 30 insertions(+), 24 deletions(-) diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index a4904f40..bb5e4c55 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -896,27 +896,6 @@ def init_bold_native_wf( ]) # fmt:skip 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([ - (echo_index, norf_source, [('echoidx', 'index')]), - (norf_source, validate_norf, [('out', 'in_file')]), - ]) # fmt:skip - - 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')]), - ]) # fmt:skip - - 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')]), - ]) # fmt:skip - dwidenoise_wf = init_bold_dwidenoise_wf( has_phase=has_phase, has_norf=has_norf, @@ -924,14 +903,41 @@ def init_bold_native_wf( ) workflow.connect([ (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 + + 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')])]) From 023c6b3f99aa7a51009182e39333eb5fb795857e Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sun, 10 Nov 2024 10:08:42 -0500 Subject: [PATCH 15/25] Write out noise images. --- fmriprep/cli/parser.py | 4 +- fmriprep/data/boilerplate.bib | 48 ++++++++++++++++++ fmriprep/interfaces/denoise.py | 13 +++-- fmriprep/workflows/bold/base.py | 1 + fmriprep/workflows/bold/denoise.py | 78 +++++++++++++++++++++++++----- fmriprep/workflows/bold/fit.py | 9 +++- fmriprep/workflows/bold/outputs.py | 43 ++++++++++++++++ 7 files changed, 178 insertions(+), 18 deletions(-) diff --git a/fmriprep/cli/parser.py b/fmriprep/cli/parser.py index 0dd9c40b..6595f0b8 100644 --- a/fmriprep/cli/parser.py +++ b/fmriprep/cli/parser.py @@ -461,8 +461,8 @@ def _slice_time_ref(value, parser): action='store', dest='thermal_denoise_method', default=None, - choices=['nordic', 'mppca'], - help='Apply NORDIC or MP-PCA denoising to the BOLD data to remove thermal noise', + 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') diff --git a/fmriprep/data/boilerplate.bib b/fmriprep/data/boilerplate.bib index 4cc802f2..d4b27d76 100644 --- a/fmriprep/data/boilerplate.bib +++ b/fmriprep/data/boilerplate.bib @@ -365,3 +365,51 @@ @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{dowdle2023evaluating, + title={Evaluating increases in sensitivity from NORDIC for diverse fMRI acquisition strategies}, + author={Dowdle, Logan T and Vizioli, Luca and Moeller, Steen and Ak{\c{c}}akaya, Mehmet and Olman, Cheryl and Ghose, Geoffrey and Yacoub, Essa and U{\u{g}}urbil, K{\^a}mil}, + journal={NeuroImage}, + volume={270}, + pages={119949}, + year={2023}, + publisher={Elsevier}, + url={https://doi.org/10.1016/j.neuroimage.2023.119949}, + doi={10.1016/j.neuroimage.2023.119949} +} + +@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/interfaces/denoise.py b/fmriprep/interfaces/denoise.py index 0865a13a..12969bb0 100644 --- a/fmriprep/interfaces/denoise.py +++ b/fmriprep/interfaces/denoise.py @@ -65,6 +65,11 @@ def _run_interface(self, 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), @@ -86,9 +91,6 @@ class DWIDenoiseInputSpec(MRTrix3BaseInputSpec): 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): @@ -115,6 +117,11 @@ class DWIDenoise(MRTrix3Base): 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' diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py index e37692f0..73088a85 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 index 606e4e76..6c7fde57 100644 --- a/fmriprep/workflows/bold/denoise.py +++ b/fmriprep/workflows/bold/denoise.py @@ -29,6 +29,7 @@ """ from nipype.interfaces import utility as niu +from nipype.interfaces.afni.utils import Calc from nipype.pipeline import engine as pe from ... import config @@ -72,18 +73,26 @@ def init_bold_dwidenoise_wf( Parameters ---------- has_phase : :obj:`bool` - True if phase data is available. False if not. + True if phase data are 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 + 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_stc_wf``) + Name of workflow (default: ``bold_dwidenoise_wf``) Inputs ------ - bold_file + 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 ------- @@ -91,13 +100,24 @@ def init_bold_dwidenoise_wf( 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__ = """\ -NORDIC or MP-PCA was applied to the BOLD data. -""" + workflow.__desc__ = ( + 'The BOLD time-series were denoised to remove thermal noise using ' + '`dwidenoise` [@tournier2019mrtrix3] ' + ) + if config.workflow.thermal_noise_estimator == 'nordic': + workflow.__desc__ += 'with the NORDIC method [@moeller2021noise;@dowdle2023evaluating].' + 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'], @@ -106,12 +126,17 @@ def init_bold_dwidenoise_wf( ) outputnode = pe.Node( niu.IdentityInterface( - fields=['mag_file', 'phase_file'], + 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( @@ -129,6 +154,7 @@ def init_bold_dwidenoise_wf( ]) # 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', @@ -152,9 +178,16 @@ def init_bold_dwidenoise_wf( 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', @@ -174,8 +207,14 @@ def init_bold_dwidenoise_wf( workflow.connect([(inputnode, complex_buffer, [('mag_file', 'bold_file')])]) # Run NORDIC + estimator = { + 'nordic': 'nordic', + 'mppca': 'Est2', + } dwidenoise = pe.Node( - DWIDenoise(), + DWIDenoise( + estimator=estimator[config.workflow.thermal_noise_estimator], + ), mem_gb=mem_gb['filesize'] * 2, name='dwidenoise', ) @@ -203,7 +242,22 @@ def init_bold_dwidenoise_wf( (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')])]) + 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 bb5e4c55..a9f31686 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -722,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.) @@ -788,7 +791,8 @@ def init_bold_native_wf( ','.join([os.path.basename(norf) for norf in norf_files]) ) config.loggers.workflow.info(norf_msg) - has_norf = bool(len(norf_files)) + # XXX: disabled until MRTrix3 implements dwi2noise + has_norf = bool(len(norf_files)) and False has_phase = False phase_files = [] @@ -862,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 @@ -907,6 +913,7 @@ def init_bold_native_wf( ('outputnode.mag_file', 'bold_file'), ('outputnode.phase_file', 'phase_file'), ]), + (dwidenoise_wf, outputnode, [('outputnode.noise_file', 'thermal_noise')]), ]) # fmt:skip if has_norf: diff --git a/fmriprep/workflows/bold/outputs.py b/fmriprep/workflows/bold/outputs.py index 70d1a4ac..533b4b7c 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( From 2679b419f527fc4b2a530e38c7951a8864cfd4aa Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sun, 10 Nov 2024 11:03:43 -0500 Subject: [PATCH 16/25] Fix variable name. --- fmriprep/workflows/bold/denoise.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fmriprep/workflows/bold/denoise.py b/fmriprep/workflows/bold/denoise.py index 6c7fde57..538217f9 100644 --- a/fmriprep/workflows/bold/denoise.py +++ b/fmriprep/workflows/bold/denoise.py @@ -110,7 +110,7 @@ def init_bold_dwidenoise_wf( 'The BOLD time-series were denoised to remove thermal noise using ' '`dwidenoise` [@tournier2019mrtrix3] ' ) - if config.workflow.thermal_noise_estimator == 'nordic': + if config.workflow.thermal_denoise_method == 'nordic': workflow.__desc__ += 'with the NORDIC method [@moeller2021noise;@dowdle2023evaluating].' else: workflow.__desc__ += ( @@ -213,7 +213,7 @@ def init_bold_dwidenoise_wf( } dwidenoise = pe.Node( DWIDenoise( - estimator=estimator[config.workflow.thermal_noise_estimator], + estimator=estimator[config.workflow.thermal_denoise_method], ), mem_gb=mem_gb['filesize'] * 2, name='dwidenoise', From 477531a60eb6b8a75c76bd9a2e2c5dc12840d749 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sun, 10 Nov 2024 11:18:33 -0500 Subject: [PATCH 17/25] Update citations. --- docs/workflows.rst | 26 ++++++++++++++++++++++++++ fmriprep/data/boilerplate.bib | 21 +++++++++++---------- fmriprep/data/tests/config.toml | 1 + fmriprep/workflows/bold/denoise.py | 2 +- 4 files changed, 39 insertions(+), 11 deletions(-) diff --git a/docs/workflows.rst b/docs/workflows.rst index 135a5ab7..e9caeee9 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/fmriprep/data/boilerplate.bib b/fmriprep/data/boilerplate.bib index d4b27d76..23280fe6 100644 --- a/fmriprep/data/boilerplate.bib +++ b/fmriprep/data/boilerplate.bib @@ -390,16 +390,17 @@ @article{tournier2019mrtrix3 doi={10.1016/j.neuroimage.2019.116137} } -@article{dowdle2023evaluating, - title={Evaluating increases in sensitivity from NORDIC for diverse fMRI acquisition strategies}, - author={Dowdle, Logan T and Vizioli, Luca and Moeller, Steen and Ak{\c{c}}akaya, Mehmet and Olman, Cheryl and Ghose, Geoffrey and Yacoub, Essa and U{\u{g}}urbil, K{\^a}mil}, - journal={NeuroImage}, - volume={270}, - pages={119949}, - year={2023}, - publisher={Elsevier}, - url={https://doi.org/10.1016/j.neuroimage.2023.119949}, - doi={10.1016/j.neuroimage.2023.119949} +@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, diff --git a/fmriprep/data/tests/config.toml b/fmriprep/data/tests/config.toml index b1b7e31b..358b7780 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/workflows/bold/denoise.py b/fmriprep/workflows/bold/denoise.py index 538217f9..d0322dad 100644 --- a/fmriprep/workflows/bold/denoise.py +++ b/fmriprep/workflows/bold/denoise.py @@ -111,7 +111,7 @@ def init_bold_dwidenoise_wf( '`dwidenoise` [@tournier2019mrtrix3] ' ) if config.workflow.thermal_denoise_method == 'nordic': - workflow.__desc__ += 'with the NORDIC method [@moeller2021noise;@dowdle2023evaluating].' + workflow.__desc__ += 'with the NORDIC method [@moeller2021noise;@vizioli2021lowering].' else: workflow.__desc__ += ( 'with the Marchenko-Pastur principal components analysis method ' From f19181eaea3ceb2b87427d00ca9562e323255275 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sun, 10 Nov 2024 11:40:29 -0500 Subject: [PATCH 18/25] Fix string. --- fmriprep/workflows/bold/denoise.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fmriprep/workflows/bold/denoise.py b/fmriprep/workflows/bold/denoise.py index d0322dad..7ed102a1 100644 --- a/fmriprep/workflows/bold/denoise.py +++ b/fmriprep/workflows/bold/denoise.py @@ -209,7 +209,7 @@ def init_bold_dwidenoise_wf( # Run NORDIC estimator = { 'nordic': 'nordic', - 'mppca': 'Est2', + 'mppca': 'Exp2', } dwidenoise = pe.Node( DWIDenoise( From f8ce52781c5897aaff11d6c407c57709f89f0f54 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Sun, 10 Nov 2024 11:58:48 -0500 Subject: [PATCH 19/25] Fix the scaling factor. --- fmriprep/workflows/bold/denoise.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fmriprep/workflows/bold/denoise.py b/fmriprep/workflows/bold/denoise.py index 7ed102a1..d2bb7450 100644 --- a/fmriprep/workflows/bold/denoise.py +++ b/fmriprep/workflows/bold/denoise.py @@ -245,7 +245,7 @@ def init_bold_dwidenoise_wf( # Apply sqrt(2) scaling factor to noise map rescale_noise = pe.Node( - Calc(expr='a*sqrt(2)', outputtype='NIFTI_GZ'), + Calc(expr='a/sqrt(2)', outputtype='NIFTI_GZ'), name='rescale_noise', ) workflow.connect([ From 2b88d701473f6d47740433fbd5deb69a0cf21aa7 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 2 Jan 2025 14:11:18 -0500 Subject: [PATCH 20/25] Try to run mppca. --- .circleci/config.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.circleci/config.yml b/.circleci/config.yml index 456fc375..4177c136 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -712,6 +712,7 @@ jobs: /tmp/data/${DATASET} /tmp/${DATASET}/fmriprep participant \ --fs-no-reconall --sloppy --write-graph \ --output-spaces MNI152NLin2009cAsym \ + --thermal-denoise-method mppca \ --mem-mb 14336 --nthreads 4 --anat-only -vv --notrack fi - run: From 7cd31e6ac2258b1c8e7de2e5853bd8e14300129e Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 2 Jan 2025 15:06:54 -0500 Subject: [PATCH 21/25] Expect noise map. --- .circleci/config.yml | 1 - .circleci/ds210_outputs.txt | 2 ++ 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/.circleci/config.yml b/.circleci/config.yml index 4177c136..456fc375 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -712,7 +712,6 @@ jobs: /tmp/data/${DATASET} /tmp/${DATASET}/fmriprep participant \ --fs-no-reconall --sloppy --write-graph \ --output-spaces MNI152NLin2009cAsym \ - --thermal-denoise-method mppca \ --mem-mb 14336 --nthreads 4 --anat-only -vv --notrack fi - run: diff --git a/.circleci/ds210_outputs.txt b/.circleci/ds210_outputs.txt index 4ee3a1e9..f18ccad2 100644 --- a/.circleci/ds210_outputs.txt +++ b/.circleci/ds210_outputs.txt @@ -41,6 +41,8 @@ sub-02/func/sub-02_task-cuedSGT_run-01_desc-coreg_boldref.json sub-02/func/sub-02_task-cuedSGT_run-01_desc-coreg_boldref.nii.gz sub-02/func/sub-02_task-cuedSGT_run-01_desc-hmc_boldref.json sub-02/func/sub-02_task-cuedSGT_run-01_desc-hmc_boldref.nii.gz +sub-02/func/sub-02_task-cuedSGT_run-01_desc-noise_boldmap.json +sub-02/func/sub-02_task-cuedSGT_run-01_desc-noise_boldmap.nii.gz sub-02/func/sub-02_task-cuedSGT_run-01_echo-1_desc-preproc_bold.json sub-02/func/sub-02_task-cuedSGT_run-01_echo-1_desc-preproc_bold.nii.gz sub-02/func/sub-02_task-cuedSGT_run-01_echo-2_desc-preproc_bold.json From ee2dc943940720916bad5f47ff3075317dd7dc2e Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 2 Jan 2025 15:07:45 -0500 Subject: [PATCH 22/25] Expect echo-wise noise maps. --- .circleci/ds210_outputs.txt | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/.circleci/ds210_outputs.txt b/.circleci/ds210_outputs.txt index f18ccad2..ca0ec3df 100644 --- a/.circleci/ds210_outputs.txt +++ b/.circleci/ds210_outputs.txt @@ -41,14 +41,18 @@ sub-02/func/sub-02_task-cuedSGT_run-01_desc-coreg_boldref.json sub-02/func/sub-02_task-cuedSGT_run-01_desc-coreg_boldref.nii.gz sub-02/func/sub-02_task-cuedSGT_run-01_desc-hmc_boldref.json sub-02/func/sub-02_task-cuedSGT_run-01_desc-hmc_boldref.nii.gz -sub-02/func/sub-02_task-cuedSGT_run-01_desc-noise_boldmap.json -sub-02/func/sub-02_task-cuedSGT_run-01_desc-noise_boldmap.nii.gz sub-02/func/sub-02_task-cuedSGT_run-01_echo-1_desc-preproc_bold.json sub-02/func/sub-02_task-cuedSGT_run-01_echo-1_desc-preproc_bold.nii.gz sub-02/func/sub-02_task-cuedSGT_run-01_echo-2_desc-preproc_bold.json sub-02/func/sub-02_task-cuedSGT_run-01_echo-2_desc-preproc_bold.nii.gz sub-02/func/sub-02_task-cuedSGT_run-01_echo-3_desc-preproc_bold.json sub-02/func/sub-02_task-cuedSGT_run-01_echo-3_desc-preproc_bold.nii.gz +sub-02/func/sub-02_task-cuedSGT_run-01_echo-1_desc-noise_boldmap.json +sub-02/func/sub-02_task-cuedSGT_run-01_echo-1_desc-noise_boldmap.nii.gz +sub-02/func/sub-02_task-cuedSGT_run-01_echo-2_desc-noise_boldmap.json +sub-02/func/sub-02_task-cuedSGT_run-01_echo-2_desc-noise_boldmap.nii.gz +sub-02/func/sub-02_task-cuedSGT_run-01_echo-3_desc-noise_boldmap.json +sub-02/func/sub-02_task-cuedSGT_run-01_echo-3_desc-noise_boldmap.nii.gz sub-02/func/sub-02_task-cuedSGT_run-01_from-boldref_to-auto00000_mode-image_xfm.json sub-02/func/sub-02_task-cuedSGT_run-01_from-boldref_to-auto00000_mode-image_xfm.txt sub-02/func/sub-02_task-cuedSGT_run-01_from-boldref_to-T1w_mode-image_desc-coreg_xfm.json From b883024c4e559ea3084ef1e37bf19cf1fb200991 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 2 Jan 2025 15:44:33 -0500 Subject: [PATCH 23/25] Run MP-PCA on ds054 too. --- .circleci/config.yml | 1 + .circleci/ds054_fasttrack_outputs.txt | 4 ++++ .circleci/ds054_outputs.txt | 4 ++++ .circleci/ds210_fasttrack_outputs.txt | 6 ++++++ 4 files changed, 15 insertions(+) diff --git a/.circleci/config.yml b/.circleci/config.yml index 456fc375..9103fa46 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -648,6 +648,7 @@ jobs: ${FASTRACK_ARG} \ --fs-no-reconall --sloppy \ --output-spaces MNI152NLin2009cAsym:res-2 anat func \ + --thermal-denoise-method mppca \ --mem-mb 14336 --nthreads 4 -vv --debug compcor - run: *check_outputs - run: diff --git a/.circleci/ds054_fasttrack_outputs.txt b/.circleci/ds054_fasttrack_outputs.txt index c9327680..87f8cacb 100644 --- a/.circleci/ds054_fasttrack_outputs.txt +++ b/.circleci/ds054_fasttrack_outputs.txt @@ -33,6 +33,8 @@ sub-100185/func/sub-100185_task-machinegame_run-01_desc-hmc_boldref.json sub-100185/func/sub-100185_task-machinegame_run-01_desc-hmc_boldref.nii.gz sub-100185/func/sub-100185_task-machinegame_run-01_desc-preproc_bold.json sub-100185/func/sub-100185_task-machinegame_run-01_desc-preproc_bold.nii.gz +sub-100185/func/sub-100185_task-machinegame_run-01_desc-noise_boldmap.json +sub-100185/func/sub-100185_task-machinegame_run-01_desc-noise_boldmap.nii.gz sub-100185/func/sub-100185_task-machinegame_run-01_from-boldref_to-auto00000_mode-image_xfm.json sub-100185/func/sub-100185_task-machinegame_run-01_from-boldref_to-auto00000_mode-image_xfm.txt sub-100185/func/sub-100185_task-machinegame_run-01_from-boldref_to-T1w_mode-image_desc-coreg_xfm.json @@ -61,6 +63,8 @@ sub-100185/func/sub-100185_task-machinegame_run-02_desc-hmc_boldref.json sub-100185/func/sub-100185_task-machinegame_run-02_desc-hmc_boldref.nii.gz sub-100185/func/sub-100185_task-machinegame_run-02_desc-preproc_bold.json sub-100185/func/sub-100185_task-machinegame_run-02_desc-preproc_bold.nii.gz +sub-100185/func/sub-100185_task-machinegame_run-02_desc-noise_boldmap.json +sub-100185/func/sub-100185_task-machinegame_run-02_desc-noise_boldmap.nii.gz sub-100185/func/sub-100185_task-machinegame_run-02_from-boldref_to-auto00000_mode-image_xfm.json sub-100185/func/sub-100185_task-machinegame_run-02_from-boldref_to-auto00000_mode-image_xfm.txt sub-100185/func/sub-100185_task-machinegame_run-02_from-boldref_to-T1w_mode-image_desc-coreg_xfm.json diff --git a/.circleci/ds054_outputs.txt b/.circleci/ds054_outputs.txt index 1340f33b..7e4f9cd6 100644 --- a/.circleci/ds054_outputs.txt +++ b/.circleci/ds054_outputs.txt @@ -43,6 +43,8 @@ sub-100185/func/sub-100185_task-machinegame_run-01_desc-hmc_boldref.json sub-100185/func/sub-100185_task-machinegame_run-01_desc-hmc_boldref.nii.gz sub-100185/func/sub-100185_task-machinegame_run-01_desc-preproc_bold.json sub-100185/func/sub-100185_task-machinegame_run-01_desc-preproc_bold.nii.gz +sub-100185/func/sub-100185_task-machinegame_run-01_desc-noise_boldmap.json +sub-100185/func/sub-100185_task-machinegame_run-01_desc-noise_boldmap.nii.gz sub-100185/func/sub-100185_task-machinegame_run-01_from-boldref_to-auto00000_mode-image_xfm.json sub-100185/func/sub-100185_task-machinegame_run-01_from-boldref_to-auto00000_mode-image_xfm.txt sub-100185/func/sub-100185_task-machinegame_run-01_from-boldref_to-T1w_mode-image_desc-coreg_xfm.json @@ -71,6 +73,8 @@ sub-100185/func/sub-100185_task-machinegame_run-02_desc-hmc_boldref.json sub-100185/func/sub-100185_task-machinegame_run-02_desc-hmc_boldref.nii.gz sub-100185/func/sub-100185_task-machinegame_run-02_desc-preproc_bold.json sub-100185/func/sub-100185_task-machinegame_run-02_desc-preproc_bold.nii.gz +sub-100185/func/sub-100185_task-machinegame_run-02_desc-noise_boldmap.json +sub-100185/func/sub-100185_task-machinegame_run-02_desc-noise_boldmap.nii.gz sub-100185/func/sub-100185_task-machinegame_run-02_from-boldref_to-auto00000_mode-image_xfm.json sub-100185/func/sub-100185_task-machinegame_run-02_from-boldref_to-auto00000_mode-image_xfm.txt sub-100185/func/sub-100185_task-machinegame_run-02_from-boldref_to-T1w_mode-image_desc-coreg_xfm.json diff --git a/.circleci/ds210_fasttrack_outputs.txt b/.circleci/ds210_fasttrack_outputs.txt index 442283e8..53337baf 100644 --- a/.circleci/ds210_fasttrack_outputs.txt +++ b/.circleci/ds210_fasttrack_outputs.txt @@ -37,6 +37,12 @@ sub-02/func/sub-02_task-cuedSGT_run-01_echo-2_desc-preproc_bold.json sub-02/func/sub-02_task-cuedSGT_run-01_echo-2_desc-preproc_bold.nii.gz sub-02/func/sub-02_task-cuedSGT_run-01_echo-3_desc-preproc_bold.json sub-02/func/sub-02_task-cuedSGT_run-01_echo-3_desc-preproc_bold.nii.gz +sub-02/func/sub-02_task-cuedSGT_run-01_echo-1_desc-noise_boldmap.json +sub-02/func/sub-02_task-cuedSGT_run-01_echo-1_desc-noise_boldmap.nii.gz +sub-02/func/sub-02_task-cuedSGT_run-01_echo-2_desc-noise_boldmap.json +sub-02/func/sub-02_task-cuedSGT_run-01_echo-2_desc-noise_boldmap.nii.gz +sub-02/func/sub-02_task-cuedSGT_run-01_echo-3_desc-noise_boldmap.json +sub-02/func/sub-02_task-cuedSGT_run-01_echo-3_desc-noise_boldmap.nii.gz sub-02/func/sub-02_task-cuedSGT_run-01_from-boldref_to-auto00000_mode-image_xfm.json sub-02/func/sub-02_task-cuedSGT_run-01_from-boldref_to-auto00000_mode-image_xfm.txt sub-02/func/sub-02_task-cuedSGT_run-01_from-boldref_to-T1w_mode-image_desc-coreg_xfm.json From 4c3654bf870eb5004d912408baba0fafd61417f1 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 2 Jan 2025 16:32:49 -0500 Subject: [PATCH 24/25] ds210 doesn't have bold_output enabled. --- .circleci/ds054_fasttrack_outputs.txt | 2 -- .circleci/ds054_outputs.txt | 2 -- .circleci/ds210_fasttrack_outputs.txt | 6 ------ 3 files changed, 10 deletions(-) diff --git a/.circleci/ds054_fasttrack_outputs.txt b/.circleci/ds054_fasttrack_outputs.txt index 87f8cacb..9891f055 100644 --- a/.circleci/ds054_fasttrack_outputs.txt +++ b/.circleci/ds054_fasttrack_outputs.txt @@ -33,7 +33,6 @@ sub-100185/func/sub-100185_task-machinegame_run-01_desc-hmc_boldref.json sub-100185/func/sub-100185_task-machinegame_run-01_desc-hmc_boldref.nii.gz sub-100185/func/sub-100185_task-machinegame_run-01_desc-preproc_bold.json sub-100185/func/sub-100185_task-machinegame_run-01_desc-preproc_bold.nii.gz -sub-100185/func/sub-100185_task-machinegame_run-01_desc-noise_boldmap.json sub-100185/func/sub-100185_task-machinegame_run-01_desc-noise_boldmap.nii.gz sub-100185/func/sub-100185_task-machinegame_run-01_from-boldref_to-auto00000_mode-image_xfm.json sub-100185/func/sub-100185_task-machinegame_run-01_from-boldref_to-auto00000_mode-image_xfm.txt @@ -63,7 +62,6 @@ sub-100185/func/sub-100185_task-machinegame_run-02_desc-hmc_boldref.json sub-100185/func/sub-100185_task-machinegame_run-02_desc-hmc_boldref.nii.gz sub-100185/func/sub-100185_task-machinegame_run-02_desc-preproc_bold.json sub-100185/func/sub-100185_task-machinegame_run-02_desc-preproc_bold.nii.gz -sub-100185/func/sub-100185_task-machinegame_run-02_desc-noise_boldmap.json sub-100185/func/sub-100185_task-machinegame_run-02_desc-noise_boldmap.nii.gz sub-100185/func/sub-100185_task-machinegame_run-02_from-boldref_to-auto00000_mode-image_xfm.json sub-100185/func/sub-100185_task-machinegame_run-02_from-boldref_to-auto00000_mode-image_xfm.txt diff --git a/.circleci/ds054_outputs.txt b/.circleci/ds054_outputs.txt index 7e4f9cd6..65adec55 100644 --- a/.circleci/ds054_outputs.txt +++ b/.circleci/ds054_outputs.txt @@ -43,7 +43,6 @@ sub-100185/func/sub-100185_task-machinegame_run-01_desc-hmc_boldref.json sub-100185/func/sub-100185_task-machinegame_run-01_desc-hmc_boldref.nii.gz sub-100185/func/sub-100185_task-machinegame_run-01_desc-preproc_bold.json sub-100185/func/sub-100185_task-machinegame_run-01_desc-preproc_bold.nii.gz -sub-100185/func/sub-100185_task-machinegame_run-01_desc-noise_boldmap.json sub-100185/func/sub-100185_task-machinegame_run-01_desc-noise_boldmap.nii.gz sub-100185/func/sub-100185_task-machinegame_run-01_from-boldref_to-auto00000_mode-image_xfm.json sub-100185/func/sub-100185_task-machinegame_run-01_from-boldref_to-auto00000_mode-image_xfm.txt @@ -73,7 +72,6 @@ sub-100185/func/sub-100185_task-machinegame_run-02_desc-hmc_boldref.json sub-100185/func/sub-100185_task-machinegame_run-02_desc-hmc_boldref.nii.gz sub-100185/func/sub-100185_task-machinegame_run-02_desc-preproc_bold.json sub-100185/func/sub-100185_task-machinegame_run-02_desc-preproc_bold.nii.gz -sub-100185/func/sub-100185_task-machinegame_run-02_desc-noise_boldmap.json sub-100185/func/sub-100185_task-machinegame_run-02_desc-noise_boldmap.nii.gz sub-100185/func/sub-100185_task-machinegame_run-02_from-boldref_to-auto00000_mode-image_xfm.json sub-100185/func/sub-100185_task-machinegame_run-02_from-boldref_to-auto00000_mode-image_xfm.txt diff --git a/.circleci/ds210_fasttrack_outputs.txt b/.circleci/ds210_fasttrack_outputs.txt index 53337baf..442283e8 100644 --- a/.circleci/ds210_fasttrack_outputs.txt +++ b/.circleci/ds210_fasttrack_outputs.txt @@ -37,12 +37,6 @@ sub-02/func/sub-02_task-cuedSGT_run-01_echo-2_desc-preproc_bold.json sub-02/func/sub-02_task-cuedSGT_run-01_echo-2_desc-preproc_bold.nii.gz sub-02/func/sub-02_task-cuedSGT_run-01_echo-3_desc-preproc_bold.json sub-02/func/sub-02_task-cuedSGT_run-01_echo-3_desc-preproc_bold.nii.gz -sub-02/func/sub-02_task-cuedSGT_run-01_echo-1_desc-noise_boldmap.json -sub-02/func/sub-02_task-cuedSGT_run-01_echo-1_desc-noise_boldmap.nii.gz -sub-02/func/sub-02_task-cuedSGT_run-01_echo-2_desc-noise_boldmap.json -sub-02/func/sub-02_task-cuedSGT_run-01_echo-2_desc-noise_boldmap.nii.gz -sub-02/func/sub-02_task-cuedSGT_run-01_echo-3_desc-noise_boldmap.json -sub-02/func/sub-02_task-cuedSGT_run-01_echo-3_desc-noise_boldmap.nii.gz sub-02/func/sub-02_task-cuedSGT_run-01_from-boldref_to-auto00000_mode-image_xfm.json sub-02/func/sub-02_task-cuedSGT_run-01_from-boldref_to-auto00000_mode-image_xfm.txt sub-02/func/sub-02_task-cuedSGT_run-01_from-boldref_to-T1w_mode-image_desc-coreg_xfm.json From 88fa7ffe812ee84863ce84976c3945ddf02ac4b2 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 3 Jan 2025 10:05:53 -0500 Subject: [PATCH 25/25] Rename variables. --- .circleci/ds210_outputs.txt | 6 --- fmriprep/interfaces/denoise.py | 42 +++++++++---------- fmriprep/workflows/bold/denoise.py | 66 +++++++++++++++--------------- fmriprep/workflows/bold/fit.py | 16 ++++---- 4 files changed, 63 insertions(+), 67 deletions(-) diff --git a/.circleci/ds210_outputs.txt b/.circleci/ds210_outputs.txt index ca0ec3df..4ee3a1e9 100644 --- a/.circleci/ds210_outputs.txt +++ b/.circleci/ds210_outputs.txt @@ -47,12 +47,6 @@ sub-02/func/sub-02_task-cuedSGT_run-01_echo-2_desc-preproc_bold.json sub-02/func/sub-02_task-cuedSGT_run-01_echo-2_desc-preproc_bold.nii.gz sub-02/func/sub-02_task-cuedSGT_run-01_echo-3_desc-preproc_bold.json sub-02/func/sub-02_task-cuedSGT_run-01_echo-3_desc-preproc_bold.nii.gz -sub-02/func/sub-02_task-cuedSGT_run-01_echo-1_desc-noise_boldmap.json -sub-02/func/sub-02_task-cuedSGT_run-01_echo-1_desc-noise_boldmap.nii.gz -sub-02/func/sub-02_task-cuedSGT_run-01_echo-2_desc-noise_boldmap.json -sub-02/func/sub-02_task-cuedSGT_run-01_echo-2_desc-noise_boldmap.nii.gz -sub-02/func/sub-02_task-cuedSGT_run-01_echo-3_desc-noise_boldmap.json -sub-02/func/sub-02_task-cuedSGT_run-01_echo-3_desc-noise_boldmap.nii.gz sub-02/func/sub-02_task-cuedSGT_run-01_from-boldref_to-auto00000_mode-image_xfm.json sub-02/func/sub-02_task-cuedSGT_run-01_from-boldref_to-auto00000_mode-image_xfm.txt sub-02/func/sub-02_task-cuedSGT_run-01_from-boldref_to-T1w_mode-image_desc-coreg_xfm.json diff --git a/fmriprep/interfaces/denoise.py b/fmriprep/interfaces/denoise.py index 12969bb0..7d303346 100644 --- a/fmriprep/interfaces/denoise.py +++ b/fmriprep/interfaces/denoise.py @@ -16,12 +16,12 @@ class _ValidateComplexInputSpec(BaseInterfaceInputSpec): - mag_file = File( + magnitude = File( exists=True, mandatory=True, desc='Magnitude BOLD EPI', ) - phase_file = File( + phase = File( exists=True, mandatory=False, desc='Phase BOLD EPI', @@ -29,8 +29,8 @@ class _ValidateComplexInputSpec(BaseInterfaceInputSpec): class _ValidateComplexOutputSpec(TraitedSpec): - mag_file = File(exists=True, desc='Validated magnitude file') - phase_file = File(exists=True, desc='Validated phase file') + magnitude = File(exists=True, desc='Validated magnitude file') + phase = File(exists=True, desc='Validated phase file') class ValidateComplex(SimpleInterface): @@ -42,12 +42,12 @@ class ValidateComplex(SimpleInterface): def _run_interface(self, runtime): import nibabel as nb - if not isdefined(self.inputs.phase_file): - self._results['mag_file'] = self.inputs.mag_file + if not isdefined(self.inputs.phase): + self._results['magnitude'] = self.inputs.magnitude return runtime - mag_img = nb.load(self.inputs.mag_file) - phase_img = nb.load(self.inputs.phase_file) + mag_img = nb.load(self.inputs.magnitude) + phase_img = nb.load(self.inputs.phase) n_mag_vols = mag_img.shape[3] n_phase_vols = phase_img.shape[3] @@ -57,8 +57,8 @@ def _run_interface(self, runtime): 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['magnitude'] = self.inputs.magnitude + self._results['phase'] = self.inputs.phase return runtime @@ -139,11 +139,11 @@ def _get_plotting_images(self): 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( + magnitude = traits.File(exists=True, mandatory=True, position=0, argstr='%s') + phase = traits.File(exists=True, mandatory=True, position=1, argstr='%s') + complex = traits.File( exists=False, - name_source='mag_file', + name_source='magnitude', name_template='%s_complex.nii.gz', keep_extension=False, position=-1, @@ -152,7 +152,7 @@ class _PolarToComplexInputSpec(MRTrix3BaseInputSpec): class _PolarToComplexOutputSpec(TraitedSpec): - out_file = File(exists=True) + complex = File(exists=True) class PolarToComplex(MRTrix3Base): @@ -251,7 +251,7 @@ class _PhaseToRadInputSpec(BaseInterfaceInputSpec): """ - phase_file = File(exists=True, mandatory=True) + phase = File(exists=True, mandatory=True) class _PhaseToRadOutputSpec(TraitedSpec): @@ -291,7 +291,7 @@ class _PhaseToRadOutputSpec(TraitedSpec): """ - phase_file = File(exists=True) + phase = File(exists=True) class PhaseToRad(SimpleInterface): @@ -338,7 +338,7 @@ class PhaseToRad(SimpleInterface): output_spec = _PhaseToRadOutputSpec def _run_interface(self, runtime): - im = nb.load(self.inputs.phase_file) + im = nb.load(self.inputs.phase) data = im.get_fdata(caching='unchanged') # Read as float64 for safety hdr = im.header.copy() @@ -352,15 +352,15 @@ def _run_interface(self, runtime): hdr.set_xyzt_units('mm') # Set the output file name - self._results['phase_file'] = fname_presuffix( - self.inputs.phase_file, + self._results['phase'] = fname_presuffix( + self.inputs.phase, 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']) + nb.Nifti1Image(data, None, hdr).to_filename(self._results['phase']) return runtime diff --git a/fmriprep/workflows/bold/denoise.py b/fmriprep/workflows/bold/denoise.py index d2bb7450..61c149b0 100644 --- a/fmriprep/workflows/bold/denoise.py +++ b/fmriprep/workflows/bold/denoise.py @@ -63,7 +63,7 @@ def init_bold_dwidenoise_wf( :graph2use: orig :simple_form: yes - from fmriprep.workflows.bold import init_bold_dwidenoise_wf + from fmriprep.workflows.bold.denoise import init_bold_dwidenoise_wf wf = init_bold_dwidenoise_wf( has_phase=True, has_norf=True, @@ -75,7 +75,7 @@ def init_bold_dwidenoise_wf( 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. + True if no-excitation (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 @@ -85,22 +85,22 @@ def init_bold_dwidenoise_wf( Inputs ------ - mag_file + magnitude BOLD series NIfTI file - phase_file + phase Phase series NIfTI file - norf_file + magnitude_norf Noise map NIfTI file - phase_norf_file + phase_norf Phase noise map NIfTI file Outputs ------- - mag_file + magnitude Denoised BOLD series NIfTI file - phase_file + phase Denoised phase series NIfTI file - noise_file + noise Noise map NIfTI file """ from niworkflows.engine.workflows import LiterateWorkflow as Workflow @@ -120,16 +120,16 @@ def init_bold_dwidenoise_wf( inputnode = pe.Node( niu.IdentityInterface( - fields=['mag_file', 'norf_file', 'phase_file', 'phase_norf_file'], + fields=['magnitude', 'magnitude_norf', 'phase', 'phase_norf'], ), name='inputnode', ) outputnode = pe.Node( niu.IdentityInterface( fields=[ - 'mag_file', - 'phase_file', - 'noise_file', + 'magnitude', + 'phase', + 'noise', ], ), name='outputnode', @@ -138,7 +138,9 @@ def init_bold_dwidenoise_wf( 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 + # TODO: Calculate the noise level from the noise volumes + # XXX: In NORDIC, noise level is estimated after scaling, phase filtering, + # and g-factor correction. noise_estimate = pe.Node( NoiseEstimate(), name='noise_estimate', @@ -148,8 +150,8 @@ def init_bold_dwidenoise_wf( 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'), + ('magnitude_norf', 'magnitude'), + ('phase_norf', 'phase'), ]), ]) # fmt:skip @@ -160,7 +162,7 @@ def init_bold_dwidenoise_wf( name='phase_to_radians_norf', ) workflow.connect([ - (validate_complex_norf, phase_to_radians_norf, [('phase_file', 'phase_file')]), + (validate_complex_norf, phase_to_radians_norf, [('phase', 'phase')]), ]) # fmt:skip combine_complex_norf = pe.Node( @@ -168,13 +170,13 @@ def init_bold_dwidenoise_wf( 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')]), + (validate_complex_norf, combine_complex_norf, [('magnitude', 'magnitude')]), + (phase_to_radians_norf, combine_complex_norf, [('phase', 'phase')]), (combine_complex_norf, noise_estimate, [('out_file', 'in_file')]), ]) # fmt:skip else: - workflow.connect([(inputnode, noise_estimate, [('mag_norf_file', 'in_file')])]) + workflow.connect([(inputnode, noise_estimate, [('magnitude_norf', 'in_file')])]) complex_buffer = pe.Node(niu.IdentityInterface(fields=['bold_file']), name='complex_buffer') if has_phase: @@ -192,19 +194,19 @@ def init_bold_dwidenoise_wf( PhaseToRad(), name='phase_to_radians', ) - workflow.connect([(validate_complex, phase_to_radians, [('phase_file', 'phase_file')])]) + workflow.connect([(validate_complex, phase_to_radians, [('phase', 'phase')])]) 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')]), + (validate_complex, combine_complex, [('magnitude', 'magnitude')]), + (phase_to_radians, combine_complex, [('phase', 'phase')]), + (combine_complex, complex_buffer, [('complex', 'bold_file')]), ]) # fmt:skip else: - workflow.connect([(inputnode, complex_buffer, [('mag_file', 'bold_file')])]) + workflow.connect([(inputnode, complex_buffer, [('magnitude', 'bold_file')])]) # Run NORDIC estimator = { @@ -230,8 +232,8 @@ def init_bold_dwidenoise_wf( name='split_complex', ) workflow.connect([ - (dwidenoise, split_magnitude, [('out_file', 'complex_file')]), - (split_magnitude, outputnode, [('out_file', 'mag_file')]), + (dwidenoise, split_magnitude, [('out_file', 'complex')]), + (split_magnitude, outputnode, [('out_file', 'magnitude')]), ]) # fmt:skip split_phase = pe.Node( @@ -239,8 +241,8 @@ def init_bold_dwidenoise_wf( name='split_phase', ) workflow.connect([ - (dwidenoise, split_phase, [('out_file', 'complex_file')]), - (split_phase, outputnode, [('out_file', 'phase_file')]), + (dwidenoise, split_phase, [('out_file', 'complex')]), + (split_phase, outputnode, [('out_file', 'phase')]), ]) # fmt:skip # Apply sqrt(2) scaling factor to noise map @@ -250,13 +252,13 @@ def init_bold_dwidenoise_wf( ) workflow.connect([ (noise_estimate, rescale_noise, [('noise_map', 'in_file_a')]), - (rescale_noise, outputnode, [('out_file', 'noise_file')]), + (rescale_noise, outputnode, [('out_file', 'noise')]), ]) # fmt:skip else: workflow.connect([ (dwidenoise, outputnode, [ - ('out_file', 'mag_file'), - ('noise_image', 'noise_file'), + ('out_file', 'magnitude'), + ('noise_image', 'noise'), ]), ]) # fmt:skip diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index a9f31686..bb7e674f 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -878,7 +878,7 @@ def init_bold_native_wf( outputnode.inputs.metadata = metadata denoisebuffer = pe.Node( - niu.IdentityInterface(fields=['bold_file', 'phase_file']), + niu.IdentityInterface(fields=['bold_file', 'phase']), name='denoisebuffer', ) boldbuffer = pe.Node( @@ -908,12 +908,12 @@ def init_bold_native_wf( mem_gb=mem_gb, ) workflow.connect([ - (validate_bold, dwidenoise_wf, [('out_file', 'inputnode.mag_file')]), + (validate_bold, dwidenoise_wf, [('out_file', 'inputnode.magnitude')]), (dwidenoise_wf, denoisebuffer, [ - ('outputnode.mag_file', 'bold_file'), - ('outputnode.phase_file', 'phase_file'), + ('outputnode.magnitude', 'bold_file'), + ('outputnode.phase', 'phase'), ]), - (dwidenoise_wf, outputnode, [('outputnode.noise_file', 'thermal_noise')]), + (dwidenoise_wf, outputnode, [('outputnode.noise', 'thermal_noise')]), ]) # fmt:skip if has_norf: @@ -922,7 +922,7 @@ def init_bold_native_wf( workflow.connect([ (echo_index, norf_source, [('echoidx', 'index')]), (norf_source, validate_norf, [('out', 'in_file')]), - (validate_norf, dwidenoise_wf, [('out_file', 'inputnode.norf_file')]), + (validate_norf, dwidenoise_wf, [('out_file', 'inputnode.magnitude_norf')]), ]) # fmt:skip if has_phase: @@ -931,7 +931,7 @@ def init_bold_native_wf( workflow.connect([ (echo_index, phase_source, [('echoidx', 'index')]), (phase_source, validate_phase, [('out', 'in_file')]), - (validate_phase, dwidenoise_wf, [('out_file', 'inputnode.phase_file')]), + (validate_phase, dwidenoise_wf, [('out_file', 'inputnode.phase')]), ]) # fmt:skip if has_phase and has_norf: @@ -943,7 +943,7 @@ def init_bold_native_wf( 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')]), + (validate_phase_norf, dwidenoise_wf, [('out_file', 'inputnode.phase_norf')]), ]) # fmt:skip else: workflow.connect([(validate_bold, denoisebuffer, [('out_file', 'bold_file')])])