Write out noise images.

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')

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}
+}

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
     <https://mrtrix.readthedocs.io/en/latest/reference/commands/dwidenoise.html>
 
+    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'

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'),
             ]),

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,32 +73,51 @@ 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
     -------
     mag_file
         Denoised BOLD series NIfTI file
     phase_file
         Denoised phase series NIfTI file
+    noise_file
+        Noise map NIfTI file
     """
     from niworkflows.engine.workflows import LiterateWorkflow as Workflow
 
     workflow = Workflow(name=name)
-    workflow.__desc__ = """\
-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

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:

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(