Dear MRtrix team,
I am currently using dwidenoise and I have a question about the exact meaning of the output -noise.

My guess is that it represents the estimated spread parameter of the underlying noise distribution. For instance: consider Rician-distributed magnitude measurements M obtained as

M = sqrt( ( A + Nreal )^2 + Nimag^2 )

where A is the underlying signal, while Nreal and Nimag are the noise in the real and imaginary channels (zero mean and standard deviation S). Then the output of option -noise is an estimate of S.

You are right, the option -noise indeed outputs a map of the estimated noise level across the image.

However, keep in mind that MP-PCA denoising operates under the assumption of Gaussian noise. Therefore, the output -noise will be an estimate of the standard deviation S of additive Gaussian noise in the magnitude images. For sufficient SNR, this is a good approximation of the Rician noise distribution. You can then correct for the remaining Rician bias with the method of moments. All the details are discussed in this paper.

In continuing with this thread, I have a follow-up question.

I have applied dpi-denoising script on a DWI with 35 gradient directions, however the noise map is only including one volume. I was expecting the same number of noise maps as the input. Would you please correct me if I’m wrong?

You did an excellent job posting your question to the forum and since most of us also have email notifications activated there’s no need to tag us personally.

About your question: as discussed earlier in this thread the option -noise will output the estimated noise level. This is indeed a single volume, as MP-PCA is directly predicated on the assumption of homoscedastic noise (constant noise level) within a local patch and across all DW volumes.

The noise level map is not the be confused with the residuals between the input and the denoised output. Indeed, the main output of dwidenoise is the 4-D denoised image with exactly the same dimensions as the input. The residuals can be calculated with mrcalc.

While people already discussed it some time ago, seems like the closest thread to ask. Anyway, as mentioned previously, one needs to fix the bias afterwards, but doing so requires a calculation based on SNR. From personal experience on phantoms, using the mean signal value from a single shell works decently well to compute the correction factor.

However, for multishell data, that won’t really work anymore. The original mrm paper mention correcting directly the signal from the kept eigen value, however that does not really apply here as one only gets a 3D map out of it. So, how do people go about correcting the maps usually in case of multishell data?

Yes indeed, Rician bias correction is not yet implemented in MRtrix. It’s on my todo list, but that list is only getting longer… If anyone wants to give it a crack (?), I’m happy to support.

I’m not entirely sure I understand your question, but I think the confusion may relate to this discussion. dwidenoise estimates the noise level, which is constant for all b-values and gradient directions. SNR is only affected in different shells because the signal attenuates, not the noise level.

The last part is exactly the issue, the SNR (or at least the signal part of it) drops with the bvalue/diffusivity. The bias correction part work as a fixed point formula of SNR, hence to initialize it you need to give it a value for the magnitude SNR. Problem is this is a voxelwise method, but we have a 4D volume of course, making it ill defined as to what is the base magnitude SNR to start the algorithm from.

Just wondering what libs are you guys using for computing special functions? I found that depending on what low level fortran/old C algorithms are bundled some are just not precise enough for the generalised version of the correction factor.

My understanding of the MRM2016 paper is that Koay’s inversion was applied on all data points independently, i.e., SNR would be initialised with the denoised signal (4-D) divided by the noise map (3-D). Alternatively, nothing is stopping you from initialising with the SNR per shell either, if you process the data per shell with the initial SNR set to the mean signal in the shell divided by the (global) noise level.

All our code uses STL and Eigen. So far we have not run into precision issues with special functions.

Well, that’s the thing, we may be talking about two different thing: the original koay paper mentions a few equations, which can be used to correct the bias for both the signal mena and the noise std using as a reference the magnitude SNR for each voxel. To do it properly as in the mppca mrm from 2015, you need, once you have split the eigen values, to use the lower half as the 3D noise map and apply it for correcting the signal on the rest of the denoised 4D maps. Doing so, you end u with a voxelwise 4D correction factor, which you can apply on both the 4D data and the 3D noise maps, although it would give you a stack of 4D noise maps after correction.

I guess I can live with that, although then in the end it would be more sensible from the physics point of view to have a single 3D noise map in the end, since it would mean each voxel through time is subject to the same contamination effet, instead of reducing back the 4D map with say the median value (although I would expect it to be pretty close in most cases), so maybe it’s a good enough workaround.

Can’t find the bessel function (they have the spherical version in eigen though, not sure how it applies versus the regular one) in those special functions’ libs, but anyway, I’ll stay with my last week quick R core wrapper, it’s overkill for now but it gets stuff done, thanks for the tip.