search

MRtrix3 / mrtrix3

MRtrix3 provides a set of tools to perform various advanced diffusion MRI analyses, including constrained spherical deconvolution (CSD), probabilistic tractography, track-density imaging, and apparent fibre density
http://www.mrtrix.org
Mozilla Public License 2.0
293 stars 180 forks source link

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

Open tsalo opened 1 day ago

tsalo commented 1 day ago

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 commented 1 day ago

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:

tsalo commented 23 hours ago

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 commented 16 hours ago

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.