ENH: Add MP-PCA/NORDIC denoising through dwidenoise

[tsalo]

Closes #3308.

Potential blockers

• dwidenoise does not support NORDIC yet, but there are plans to add it (see ENH dwidenoise: Support NORDIC method and/or allow noise scans MRtrix3/mrtrix3#3031). In the meantime, MP-PCA is the only option.
• dwidenoise doesn't accept a precalculated noise map yet, which limits reproducibility if users want to apply the denoising step using the output noise map using precalculated data.
• MRTrix3 doesn't have a dwi2noise command to estimate the noise map from no-RF volumes, so I've disabled no-RF processing for now.
• Tracking lost degrees of freedom is important and dwidenoise has a -rank option to output a DOF map, but that feature hasn't been released yet. See ENH dwidenoise: Optionally output degrees of freedom map MRtrix3/mrtrix3#2312.
• I'm not sure yet what figures we should generate for the denoising procedure. I think stat-maps of the noise map, DOF map, and possibly variance-explained map (which we'll need to generate) would be good. Maybe a histogram of the noise values?
• The MRTrix3 team recommends concatenating data across echoes before running dwidenoise, and they ultimately plan to support 5D NIfTIs, but my implementation currently just applies the denoising step independently to each echo. See ENH dwidenoise: Support >4D data MRtrix3/mrtrix3#3021.
• I need to add an appropriate filename for the noise map to the niworkflows specification (pending Add boldmap suffix to filename patterns niworkflows#899).

Changes proposed in this pull request

• Create a workflow to apply MP-PCA or NORDIC denoising to BOLD data. This is done right before slice timing correction in init_bold_native_wf.
• Add a new parameter, --thermal-denoise-method to specify which approach to use, if any.
• Add "phase" and "norf" options to the ignore parameter to let users skip those elements of the thermal denoising procedure.

[codecov]

Codecov Report

Attention: Patch coverage is 53.61702% with 109 lines in your changes missing coverage. Please review.

Project coverage is 71.07%. Comparing base (08fecee) to head (f8ce527).

Files with missing lines	Patch %	Lines
fmriprep/workflows/bold/fit.py	35.71%	24 Missing and 12 partials ⚠️
fmriprep/interfaces/denoise.py	70.00%	33 Missing ⚠️
fmriprep/workflows/bold/denoise.py	40.81%	24 Missing and 5 partials ⚠️
fmriprep/workflows/bold/outputs.py	0.00%	6 Missing and 1 partial ⚠️
fmriprep/utils/bids.py	60.00%	2 Missing and 2 partials ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #3395      +/-   ##
==========================================
- Coverage   72.06%   71.07%   -0.99%     
==========================================
  Files          57       59       +2     
  Lines        4231     4464     +233     
  Branches      455      482      +27     
==========================================
+ Hits         3049     3173     +124     
- Misses       1067     1156      +89     
- Partials      115      135      +20     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[Lestropie]

dwidenoise doesn't accept a precalculated noise map yet, which limits reproducibility if users want to apply the denoising step using the output noise map using precalculated data.

It is intended to address this at some stage; for tracking see MRtrix3/mrtrix3#2274 and MRtrix3/mrtrix3#3035.

MRTrix3 doesn't have a dwi2noise command to estimate the noise map from no-RF volumes, so I've disabled no-RF processing for now.

I've added a note to MRtrix3/mrtrix3#3035.

Tracking lost degrees of freedom is important ...

Agree 100%; we're invested locally in being able to modulate downstream applications appropriately based on this. It would be good to:

• Export noise and rank images to the HTML for visual assessment.
• Apply estimated transforms to the rank image and aggregate, to produce an approximate signal rank per output element (ie. needs to be propagated for each output space, including surfaces)

dwidenoise has a -rank option to output a DOF map, but that feature hasn't been released yet.

There's a few branches in parallel that have proposed changes to dwidenoise that may be on the dev development branch, or are not even there yet let alone on master. If need be I can aggregate multiple such changes onto a dedicated branch that fMRIPrep can then clone.

... possibly variance-explained map (which we'll need to generate) ...

That could be done. It would be reasonably easy to add an export option to dwidenoise rather than trying to compute it after the fact. Do you want to make a feature request?

The MRtrix3 team recommends concatenating data across echoes before running dwidenoise, and they ultimately plan to support 5D NIfTIs, but my implementation currently just applies the denoising step independently to each echo.

When I've done this myself, I've concatenated data across echoes into the fourth dimension, denoised, and then separated the volumes back out into echoes. That is I think more principled given that the noise level should be identical. It's not super difficult to do. You may however want to clamp the maximal patch size: once (# echoes x # volumes) exceeds 343 the default changes to a 9x9x9 patch and that runs very slowly. This behaviour should be better with MRtrix3/mrtrix3#2742 / MRtrix3/mrtrix3#3029, which won't have such step changes in patch size.

---

Some other points to consider:

1. I found the denoising under default usage to be imperfect for fMRI data due to the large offset from zero. This seems to be possible to mitigate by either:

   • Demeaning the data prior to PCA: dwidenoise: Add option -demean MRtrix3/mrtrix3#2363
     This however would benefit from confirmation that it does not introduce any kind of bias into the noise estimator / reported rank.
   • Utilising double-precision floating-point.
     This can already be chosen at the command-line, but obviously comes with a performance penalty.
   • Using a more stable PCA method: dwidenoise: try using bidiagonal divide and conquer SVD MRtrix3/mrtrix3#2906
     This also comes with a performance penalty.

2. The dwidenoise command has a different copyright to the default MRtrix3 one. There is also a patent owned by NYU regarding the denoising of a 4D series based on the MP distribution. There's been internal MRtrix3 discussion on whether the scope of that patent is restricted to the specific rank estimator presented in the original article, and therefore does not cover the alternative estimator in use by default as published in this work, but I don't think we reached a definitive conclusion.
   Pinging @jdtournier, @dchristiaens, @jelleveraart for full transparency / invitation of opinions.

[jelleveraart]

Thanks @Lestropie for looping me in this conversation.

• I agree that accepting precomputed noise maps is an impactful feature, not only because it might accelerate the algorithm but also improve its robustness in case of correlated noise (https://doi.org/10.1162/imag_a_00049). It would further allow denoising of individual series without the need for concatenating them.

• Similarly, I wonder whether accepted the "rank" as an input might be of interest for certain application in which preserving homoscedasticity is more critical than maxizing precision.

• Concatenating imaging series with varying echos in the fourth dimension is conceptually a proper strategy under the assumption of identical noise levels. However, as suggested by @Lestropie, it is recommended to keep the patch size relatively small if the noise level is spatially varying. The 5x5x5 is commonly used. If the individual echo series have > 125 images, then denoising the individual series might actually be preferred as it does not impose any risk of masking subtle signal features across the series. An approach could be to (a) concatenate all series to estimate the noise level; (b) use the resulting noise map as an input to denoise the individual series.

• We have previously used the demeaning step, but I will confirm that there is not risk of introducing biases or lowering the robustness of the current algorithm.