ENH dwidenoise: Support NORDIC method and/or allow noise scans

[tsalo]

Is your feature request related to a problem? Please describe.
I would really like to include a NORDIC denoising step in fMRIPrep, QSIPrep, and ASLPrep, but the official implementation of NORDIC is in MATLAB and is pretty difficult to use. Would you be willing to incorporate the NORDIC algorithm in dwidenoise?

Describe the solution you'd like
A new -algorithm or -method parameter for dwidenoise that accepts a value of either mppca or nordic would be amazing.

Describe alternatives you've considered
Could dwidenoise accept a 4D noise image (either magnitude or complex) instead of a float for -noise? This would be necessary for NORDIC, but I don't know if MP-PCA could leverage the same information.

[Lestropie]

Not sure if this was prompted by the fact that I'm currently at a Hackathon working on dwidenoise improvements: #3029.

I also recently generated a meta-Issue for listing the various improvements that would be useful in this context: #3023.

I am interested in improving the efficacy of dwidenoise for both DWI and fMRI data. @BahmanTahayori has a pending manuscript on the benefits of using MP-PCA for multi-echo fMRI data analysis. We will be advocating that this be done upstream of fMRIPrep, combined with particular options to provide to TEDANA after fMRIPrep. It would however be more elegant if PCA denoising, whether MP-PCA or something else, were to be executed by fMRIPrep itself, and I'd be quite happy if I could progress dwidenoise to the point of being that choice.

The ability to have dwidenoise take as input a pre-estimated noise map, and use that to threshold the eigenvalue spectrum rather than running the MP estimator, would both bring the usage more in line with NORDIC, and work nicely with fMRIPrep's new fit-apply model, as the final estimated nose level could be written as an intermediate derivative. This is listed as #2274 and I currently have a team member working on it.

I'm not intimately familiar with the exact distinctions between NORDIC and complex MP-PCA. I believe there's something about g-factor penalty, and I don't know if that's possible to estimate without explicit calibration data. We could at least have access to the slice timing vector to know which slices were members of the same group, but I suspect that's not adequate.

Regarding how to provide that input noise information, there's multiple sources of confusion:

• The existing -noise option is an output, not an input; but as described above, there will soon be an alternative option that takes a noise map as input
• A noise map would be 3D rather than 4D
• I'm not sure how a noise map would be complex. Is this claiming different noise levels between real and imaginary components? I would not expect to see such a thing from a quadrature receiver array. It seems from looking at the code, the eigenvalues are complex but the MP distribution is fit to their magnitudes.

[tsalo]

Not sure if this was prompted by the fact that I'm currently at a Hackathon working on dwidenoise improvements

That was more of a happy coincidence, but I did look through your recent dwidenoise issues to see if there was any overlap. I didn't think there was, but I completely missed #2274.

The existing -noise option is an output, not an input; but as described above, there will soon be an alternative option that takes a noise map as input

Ah, sorry. I didn't realize that.

A noise map would be 3D rather than 4D

If possible, it would be great to be able to provide the no-excitation volumes (organized in BIDS as files with the _noRF suffix) directly, rather than having to calculate a 3D noise map ahead of time, but I understand if you consider that out of scope. The protocol I've been seeing (and using) lately includes 3 no-excitation volumes at the end of the functional scan, so it's definitely a 4D file initially.

I'm not sure how a noise map would be complex. Is this claiming different noise levels between real and imaginary components? I would not expect to see such a thing from a quadrature receiver array. It seems from looking at the code, the eigenvalues are complex but the MP distribution is fit to their magnitudes.

I don't actually know either. NORDIC accepts complex-valued BOLD data with no-excitation volumes at the end, so they definitely have a phase component, but NORDIC could be discarding the phase component of the noise volumes within the workflow. Even if that is the case, it would be nice to be able to provide the complex-valued noise volumes instead of having to organize the BOLD data and noise volumes differently, but I understand if that's out of scope.

On the other hand, I don't know how to convert the phase part of the noise volumes to radians, which might make it hard to create the complex-valued version of the noise volumes. NORDIC just accepts the magnitude and phase data separately, without any scaling.

[Lestropie]

If possible, it would be great to be able to provide the no-excitation volumes (organized in BIDS as files with the _noRF suffix) directly, rather than having to calculate a 3D noise map ahead of time, but I understand if you consider that out of scope.

I don't think it's "out of scope", I would rather advocate for handling it differently.

If there's the capability to provide a noise level estimate as input, that noise level can be estimated through any means. My main goal in doing so is to estimate the noise level based on the MP distribution, filter that estimate (eg. smoothing, outlier detection & replacement), and then apply that noise level estimate for denoised series reconstruction. However if one has acquired data with no RF, you can do similar by explicitly estimating a noise level map from those data.

So while a higher-level workflow for data denoising may need to branch depending on what data are available / how one chooses to filter the noise map, and that may include splitting a 4D series based on the presence / absence of RF excitation, I would personally keep the dwidenoise interface very simple, just taking a noise map as input.

Is there a publicly visible dataset that includes the no-RF volumes at the end of a 4D series? I'm curious to know how this is encoded in the DICOMs and whether our software needs to do something smarter with it. I've also never seen a pulse sequence with this capability; do you know if this is GE only?

I'm not sure how a noise map would be complex.

NORDIC accepts complex-valued BOLD data with no-excitation volumes at the end, so they definitely have a phase component, but NORDIC could be discarding the phase component of the noise volumes within the workflow.

I think we're describing different things. A noise level map I expect to be real-valued. What I think you're referring to is the volumes to be used for noise mapping. Those can be (and ideally should be) complex just as can / should be the input data to be denoised; but if an explicit calculation is performed using those data as input to estimate the noise level, I expect that to yield a real-valued result.

I don't know how to convert the phase part of the noise volumes to radians, which might make it hard to create the complex-valued version of the noise volumes.

This is a common existing step for denoising of complex DWI data. Some tools will look at the minimum and maximum values in a phase image and linearly transform the intensities so that these lie in for instance the [-pi, pi] range, though there can be some small error in doing it this way if one is not careful. For Siemens data phase images are known to be encoded in the [-4096, 4094] range. dwidenoise expects data of type "complex float"; mrcalc has a -polar option that converts separate magnitude and phase images into the real & imaginary components of a native complex floating-point type.

[tsalo]

If there's the capability to provide a noise level estimate as input, that noise level can be estimated through any means.

I see, so you'd prefer to have a separate step to create the noise map from the no-RF volumes. Would you be willing to include that as a separate command in mrtrix3 or would you prefer that be handled by a different package?

Is there a publicly visible dataset that includes the no-RF volumes at the end of a 4D series?

I have a multi-echo BOLD dataset on OpenNeuro with no-RF volumes included. I know of a few other folks who are acquiring similar data (e.g., Damien Fair's group at UMinn, a number of folks at WUSTL, Cesar Caballero-Gaudes, and Dan Handwerker), but I don't think anyone else has had a chance to share their data yet.

I'm not aware of any DWI datasets with no-RF volumes, but I am pretty sure that NORDIC is designed to leverage no-RF volumes in DWI as well.

I'm curious to know how this is encoded in the DICOMs and whether our software needs to do something smarter with it. I've also never seen a pulse sequence with this capability; do you know if this is GE only?

As far as I know, this is only (or at least most visibly) implemented in the Siemens CMRR-MB sequence. I've chatted with a number of other folks and have also looked at the DICOMs myself, and the consensus seems to be that there's no info in the DICOM that can be used to distinguish the regular volumes from the no-RF volumes. I typically include a post-heudiconv curation script to split out the no-RF volumes into files with the _noRF suffix. Unfortunately, I can't share our DICOMs (sorry about that).

Perhaps the CMRR folks would be willing to modify the sequence to encode that info in the DICOMs though.

A noise level map I expect to be real-valued. What I think you're referring to is the volumes to be used for noise mapping.

Ah, that makes perfect sense.

This is a common existing step for denoising of complex DWI data.

I was concerned about rescaling the signal data and no-RF data separately in case that might result in different scaling, but I could probably concatenate them back together for the rescaling step.

[Lestropie]

Would you be willing to include that as a separate command in mrtrix3 or would you prefer that be handled by a different package?

It will be possible to estimate via the existing MP-PCA method. I am considering the prospect of having a dwi2noise command that does almost everything dwidenoise does, it just has as its compulsory output argument the estimated noise level rather than a denoised 4D series. That way one could run that command, do any kind of filtering of the resulting noise map they deem suitable, and then run dwidenoise with an input noise map.

In the case of the no-RF data, it would require a separate implementation to estimate the noise map as written in the NORDIC manuscript. It's not the most complex thing in the world, but it'd need to be written regardless. The g-factor penalty may need to be computed differently given that sensitivity profiles would not be available. So as a first pass I'd probably be leaning on the MP-PCA estimator, and use of the no-RF data via the method described in the NORDIC manuscript would be a future augmentation.

I was concerned about rescaling the signal data and no-RF data separately in case that might result in different scaling, but I could probably concatenate them back together for the rescaling step.

If you know a priori how the phase images are numerically encoded in DICOM, then the multiplier to get to radians is fixed and calculable a priori also. There is a small risk with the data-driven approach if one series contains both a voxel with value -4096 and a voxel with value +4094 whereas the other does not; but I wouldn't conceptualise that as an error of consistency so much as that for at least one of the two, the conversion is wrong. I've a suspicion that the information about how phase is encoded may be buried somewhere within private fields, in which case conversion to radians could occur at the import stage (#3009).

[tsalo]

I've opened a PR to fMRIPrep to add a denoising step with dwidenoise: nipreps/fmriprep#3395. There are a few blockers, but I'd love your input on the overall structure.

[Lestropie]

Thanks for the ping; see response at nipreps/fmriprep#3395.

[tsalo]

SInce I've translated the NORDIC MATLAB implementation to Python line by line (often without fully understanding the steps) (see here), I was thinking I could describe what it does to see what elements, if any, would need to be incorporated into dwidenoise (beyond the SV thresholding of course).

I've struck through the steps that explicitly happen outside of dwidenoise.

1. Scale the phase data to -pi to pi.
2. Combine the magnitude and phase data into complex-valued data.
3. Normalize the complex data by the minimum non-zero magnitude in the first volume.
4. Calculate "meanphase" (unused) as the mean of the complex data (ignoring noise volumes).
5. Filter the phase data.
6. Multiply the complex data by the exponential of the angle of the filtered phase to get "k-space" data.
   • Just to be clear, I don't think this is effective k-space data, but that's what the variable names in the original code seem to imply.
7. Limit the number of volumes to 90 or fewer for g-factor estimation.
8. Estimate the g-factor map from the k-space data.
   • AFAICT, the g-factor map is the noise variance map from running MP-PCA on this version of the data.
   • Discard this version of the k-space data.
9. Identify if there are any zero elements in the g-factor map.
10. Normalize the complex data by the g-factor map.
11. Recalculate the k-space data using the normalized complex data and the meanphase.
12. Estimate the noise level from the noise volumes of the k-space data.
13. If temporal_phase is 3, perform a secondary step for filtered phase with residual spikes.
14. Multiply the k-space data by the exponential of the angle of the filtered phase.
15. If there are zero elements, fill them in with random noise.
16. Create an NVR threshold by creating random arrays in the shape of the Casorati matrix, running SVD on them, and retaining the first singular value. Take the mean across 10 iterations. Scale the NVR threshold using an error factor (default is 1) and the measured noise level. Plus sqrt if data are complex.
17. Patch-based denoising. The only NORDIC-specific thing here, AFAICT, is the actual threshold, although the patch shape is different between the g-factor estimation and NORDIC steps.
18. Rescale the denoised data by the g-factor map, then the filtered phase, then the "absolute scale" from step 3.
19. For the magnitude data, if full_dynamic_range is True, scale the data using a "gain level".
20. For the phase data, scale the data by the original range and center.

[Lestropie]

I'll make a first attempt at responses; some of these might require links to code lines and further discussion, or more time put into trying to follow their manuscript.

4. Calculate "meanphase" (unused) as the mean of the complex data (ignoring noise volumes).

I believe what's happening here is something like this.

The on-scanner reconstruction has to figure out the right way to, for each slice, combine complex data across the multiple receiver elements. Depending on how this is computed by the sequence, it's possible that even in the absence of field inhomogeneity / T2* effects / sample motion, the complex signal throughout the sample may still not be strictly real within that slice. Exactly what it will look like will depend on the sequence recon and on the receive array. For some analyses this may not matter, but for denoising it could artificially inflate the signal rank, especially if the phase offset is inconsistent between adjacent slices that contribute to the same Casorati matrix. Note that it is a known limitation with the CMRR sequence pre-R017 that the net reconstructed phase can be highly inconsistent between adjacent slices.

The approach they are seemingly using here is based on the observation that while the phase map of each slice of each volume may vary, computing the mean phase map for that same slice across all acquired volumes will give a reasonable estimate of the overall "uninteresting" phase pattern resulting from the details of the sequence multi-coil reconstruction. Explicitly computing this mean, and subtracting the resulting phase angles from the data to be denoised, therefore assists in reducing signal rank.

I can see this being more effective for fMRI data than DWI data. In the former case it may also help to reduce signal rank in areas of strong susceptibility artifacts. If it were to be used for DWI data I think that using exclusively the b=0 volumes would make it more effective, given b!=0 volumes can have phase differences due to microscopic motion. However I also think that if other mechanisms for phase demodulation are effective (see below), this step may be redundant.

5. Filter the phase data.

There may be some confusion regarding what's happening here, further evidenced by some code comments. The goal isn't to 1) get phase and 2) filter it. The goal is to demodulate gross phase differences. What you ideally want for denoising is complex data where the bulk signal amplitude is along the positive real axis, you have Gaussian noise +/- of that in the real axis, and you have Gaussian noise of equivalent variance in the imaginary axis. The goal is to compute an appropriate phase map such that when the empirical data are demodulated by that map, that is the result you observe. You can't just demodulate by the empirical phase, as then you just get a zeroed imaginary component and you're effectively denoising magnitude-transformed data. So you need some kind of "low-pass-filtered" version of the phase. Note however that approaches to do this don't operate directly on phase data; they perform some manipulation of complex data, and then just the phase of that result is used to demodulate / re-apply the modulation.

My dwidenoise feature branch now has two solutions:

1. An enhanced version of the linear phase ramp estimation shown in this work. For each slice for each timepoint, the location and phase of the peak in "k-space" (N.B. Fourier space; see below) is interpreted as the in-plane direction, spatial frequency, and phase offset, of a linear phase ramp for demodulation of the empirical data in that slice. See: dwidenoise: Low-frequency phase removal #3037. I improved it slightly by doing a Newton optimisation on the l-space data rather than just selecting the maximum within an upsampled k-space.

2. A Tukey-FFT-filtered version of the empirical data. I made the choice of window shape and default width based on this work. I don't know how much investment was made in tuning that filter, or how it may behave in comparison to what NORDIC has in place.

6. • Just to be clear, I don't think this is effective k-space data, but that's what the variable names in the original code seem to imply.

I believe that anywhere "k-space" is referenced in any work that does not explicitly denoise data prior to multi-coil reconstruction, you should just think of it as "FFT-transformed data".

7. Estimate the g-factor map from the k-space data.
   AFAICT, the g-factor map is the noise variance map from running MP-PCA on this version of the data.

If that's what they're doing, then I suspect that referring to such data as strictly a g-factor map is potentially erroneous. In particular, the "Prescan normalise" filter commonly utilised on Siemens scanners produces a strong noise nonstationarity, since it's essentially a multiplicative B1- bias field correction. So talking more generally, this is something like a "noise level map".

If that interpretation is correct, then my current solution facilitates (though does not yet expedite) a similar approach. One can take a noise level map estimated through one run of MP-PCA (eg. new dwi2noise command), and use that to modify a subsequent dwidenoise execution, scaling voxel intensities to correct for within-patch noise level variation. Eventually I want to integrate into dwidenoise an iterative approach that integrally estimates and refines this map.

18. Rescale the denoised data by the g-factor map, then the filtered phase, then the "absolute scale" from step 3.

In my latest dwidenoise code at #3029 I chose to collapse this all into the concept of "preconditioning". Demeaning, phase demodulation, and noise nonstationarity, are all estimated from the data, their inverse is applied to the data prior to denoising, and then they are re-applied to the denoised data.

I'm a little confused by the choice of operation 3. though. To me the minimal non-zero value in the first volume feels arbitrary. I can imagine something like trying to get hard-coded parameters to result in comparable manifest behaviour even where different input datasets encode numerical values across different orders of magnitude, but even then it feels like the wrong operation. But I could well be missing something.