Noise estimation

[JohnAshburner]

Noise estimation with the Rice mixture sometimes uses the variance for the wrong Rice distribution - particularly if image noise is closer to non-central chi. I'd guess the code uses the one with the smaller mean parameter, but this is sometimes not the background class because a parameter of zero with a large variance might be a better fit to the intensity distribution within the head. Maybe it would be better to select the one with the smaller expectation:

% https://en.wikipedia.org/wiki/Rice_distribution t = -nu.^2./(2sig2); % $L_{1/2}(t)$ is computed by: L = (1-t).besseli(0,-t/2,1)... -t.besseli(1,-t/2,1); % Mean (expectation) of Rice distribution f = sqrt(pisig2/2)*L

[JohnAshburner]

Sorry. Bad use of markdown. This is better...

% https://en.wikipedia.org/wiki/Rice_distribution
t  = -nu.^2./(2*sig2);
% $L_{1/2}(t)$ is computed by:
L  = (1-t).*besseli(0,-t/2,1)...
        -t.*besseli(1,-t/2,1);
% Mean (expectation) of Rice distribution
f  = sqrt(pi*sig2/2)*L

[brudfors]

Hi @JohnAshburner

That should actually already be the case:

https://github.com/balbasty/nitorch/blob/b8728103cc33285509e6f157f20265f57b38dad4/nitorch/vb/mixtures.py#L481

(I just implemented it the way you are doing it in spm_noise_estimate)

[balbasty]

Maybe I am not using your function correctly Mikael. In estimate_noise, it seems to me that if mu_noise is not provided, the class with smallest variance is assumed to be the noise class. What should I provide as "mu_noise"? John suggests to use the class with smallest mean instead of smallest variance.

[brudfors]

Sorry, I reopened.

I think I am a bit confused. JA wrote: "I'd guess the code uses the one with the smaller mean parameter", but I guess he meant variance, as that currently defines the noise class (the one with the smallest variance).

It also provides the mean - calculated the way JA suggests above - so we should base setting the noise class on the mean then?

[balbasty]

We currently use the class with smallest sd (I think). Alternatives are: using the class with smallest mean parameter (i.e. the parameter returned by the RMM class), or the class with smallest expected mean. John suggests using the latter because when the data is in reality chi distributed instead of rice, the order of the mean parameters can be different from the order of the expected value (?)

[balbasty]

Oh I misread you. So currently the RMM class returns the expected value of the distribution, not the mean parameter?

[brudfors]

Yes, it should return exactly what JA proposes.

[balbasty]

Ok I'll see if the switch works and will maybe rename the output of estimate_noise so it's clearer. Thanks for the help!

[JohnAshburner]

I haven't actually checked the code, but for some reason Klara was getting very high noise estimates compared on a dataset compared to what we were getting with a mixture of two chi distributions. My guess was that it was choosing the Rice distribution with the smaller mean parameter. It seems there must be another reason for the difference.

[JohnAshburner]

Taking the smaller variance should also be fine. For the t1w_mfc_3dflash_v1i_R4_0015/anon_s2018-02-28_18-26-190921-00001-00224-1.nii file from the MPM example dataset, I get a variance estimate of 320.3 using spm_noise_estimate.m (Rice mixture) and 539.0 from fitting a mixture of chi distributions. I'm not sure what nitorch would give.

Klara's nitorch experiments gave variances of 2856.5, 3392.1 and 2847.7 for the three runs, whereas our chi mixture model implementation gave an overall variance estimate of 386.8. This makes quite a difference to the strength of the denoising.

[brudfors]

Where can I find the example dataset? I will try it myself :)

[JohnAshburner]

The link from https://hmri-group.github.io/hMRI-toolbox/ takes you to https://hmri-group.github.io/hMRI-toolbox/.

[JohnAshburner]

Meant to say it takes you to https://owncloud.gwdg.de/index.php/s/iv2TOQwGy4FGDDZ.

[brudfors]

nitorch's noise estimate gives a variance of 322.9 for t1w_mfc_3dflash_v1i_R4_0015/anon_s2018-02-28_18-26-190921-00001-00224-1.nii, so very simillar to SPM I'd say.

[brudfors]

Do you know what images Clara used for her experiments?

[JohnAshburner]

All very strange. She used the example hMRI dataset. Perhaps she forked an old version of nitorch, or a bug has since been introduced somewhere. My own Python skills are not up to the job of tracking down the problem.

I think I'm probably wasting your time, so maybe the issue should be closed.

[brudfors]

There might still be a bug, maybe we can add Clara to this discussion, and I could ask her? Do you know if she has a github username?

[brudfors]

I will try to compare on the hMRI dataset!

[JohnAshburner]

See https://github.com/plutoniusz

[plutoniusz]

nitorch's noise estimate gives a variance of 322.9 for t1w_mfc_3dflash_v1i_R4_0015/anon_s2018-02-28_18-26-190921-00001-00224-1.nii, so very simillar to SPM I'd say.

Is there a proper way to print the noise estimate values? I currently just added a "print (sd_noise, sd_not_noise, mu_noise, mu_not_noise)" line in the estimate_noise and for the dataset https://owncloud.gwdg.de/index.php/s/iv2TOQwGy4FGDDZ/download?path=%2F&files=hmri_sample_dataset_with_maps.zip I get values around [50, 300, 70, 600] (roughly). So sd_noise^2 is about order higher than expected. Maybe I am just looking at the wrong thing?

[brudfors]

Hi @plutoniusz

Thanks for looking into this :)

This is how I do it, and the results I get:

matlab

p = 'anon_s2018-02-28_18-26-190921-00001-00224-1.nii';

sd_noise = spm_noise_estimate(p);

disp(sd_noise)  % 17.8982

python

from nitorch.tools.img_statistics import estimate_noise

p = "anon_s2018-02-28_18-26-190921-00001-00224-1.nii";

sd_noise = estimate_noise(p)[0]

print(sd_noise)  # 17.9737

As you can see, for that image, the results are quite simillar. I have not tested it for all images in the sample dataset though.

[plutoniusz]

Hi,

thank you. Yes, that works for me for the whole dataset. I have now identified that the values change once I also use nitorch.tools.qmri.io.GradientEchoMulti.from_fnames() that I used before creating maps with nitorch.tools.qmri.relax.greeq(). Based on what you sent I added that step:

from nitorch.tools.img_statistics import estimate_noise
import os
import nitorch.tools.qmri.io as qio

cwd = os.getcwd()
p = os.path.join(cwd, "MPM/t1w_mfc_3dflash_v1i_R4_0015")
#t1w = str(os.path.join(p, 'anon_s2018-02-28_18-26-190921-00001-00224-1.nii'))

t1w = []
for filename in os.listdir(str(p)):
    if filename.endswith(".nii"):
        t1w.append(str(os.path.join(p, filename)))
    else:
        continue

t1w = qio.GradientEchoMulti.from_fnames(t1w)

for e in t1w:
    print(e)
    sd_noise = estimate_noise(e.fdata())[0]
    print(sd_noise)  
# GradientEcho(shape=(280, 320, 224)) 
# tensor(59.0657, dtype=torch.float64)

Is that expected? I think that the parameter maps created with the qmri.greeq() looked more denoised than expected although I am not sure.

[brudfors]

The qmri tooling is not something that I have used. @balbasty do you perhaps know what might be going on here?

[balbasty]

The one big difference is that GradientEchoMulti uses rand=True by default when loading the data (so adds uniform noise scaled by the nifti slope tag). @plutoniusz can you try the same snippet but with e.fdata(rand=False)?

[plutoniusz]

Ok, yes with e.fdata(rand=False) the values are back to ~18. In nitorch.tools.qmri.relax._mpm._preproc.preproc() the rand=True:

dat = echo.fdata(**backend, rand=True, cache=False)
sd0, sd1, mu0, mu1 = estimate_noise(dat)

So if I estimate the noise outside of nitorch I should take into account that GradeintEchoMulti adds noise, correct? I am not sure why and where it happens?

[balbasty]

I don't know what's best to do. The point is that when the data is stored in int, there is some uncertainty about the actual value of the original data: the integer value 16 can originate from any floating point value in [16, 17).

In @JohnAshburner's code, the histogram's bin width takes the slope tag into account, so that it does not matter:

• https://github.com/spm/spm12/blob/master/toolbox/Longitudinal/spm_noise_estimate.m

In @brudfors's code, the bin width is (max - min)/1024, whether the input data is int or float.

• https://github.com/balbasty/nitorch/blob/master/nitorch/tools/img_statistics.py#L87

I didn't expect that adding the scaled uniform noise would make such an impact... What I could do is make noise_estimate also accept MappedArray inputs and, in that case, take the slope into account when computing the histogram.

[balbasty]

Update! Both John and Mikael exclude voxels with value zero from the histogram computation. When I add uniform noise, all zeros become nonzero and are included in the histogram. In the hMRI demo dataset, that makes a huge difference since images are defaced and therefore have lots of artificial zeros. Maybe I should only add noise to nonzero voxels, to keep the implicit "missing data" mask (at least in the qMRI bits where we know we're working with MRIs, so "zero" means something). I'll modify the qMRI io stuff so that this is dealt with in a nicer way and is transparent for you.

@plutoniusz @brudfors

[plutoniusz]

@balbasty thank you for the explanation

[balbasty]

This should be fixed by this commit: 5904818 (and 38030fe, forgot to scale the noise) Basically, the GradientEcho class now treats zeros as missing values and converts them to NaNs. Since we were already masking NaNs out in the fitting code, the rest should work as before.

I've also made estimate_noise use the slope tag to define the bin width when the input is a file (same as in John's code).

[brudfors]

That sounds great @balbasty :) Let us know how it goes @plutoniusz, and if it all works for you we can close this issue.

[plutoniusz]

I think it works, the values are around 18 for this dataset

[balbasty]

Thanks for the comments on the code Klara! I was messy, should be fixed now.

[JohnAshburner]

I was sure I had the following working (line 128 of nitorch/tools/img_statistics.py): dat = dat.fdata(rand=True, missing=0)

An earlier workaround simply involved inserting W[0]=0. at about line 152, but this is a bit messy.

[balbasty]

I think that the missing keyword only exists in the qmri-specific io, but not in the general io module. Maybe I should also add this missing stuff in the main io?

[JohnAshburner]

Thanks for re-opening. That would be helpful Yael. I got a bit confused when I found the missing option didn't work in some other tests I was trying - although I have to admit I didn't spend too long trying to figure it all out.

[balbasty]

What do you think we should do with values considered as missing (apart from not sampling noise): keep their original value / replace them with 0 / replace them with NaN / ... ? I feel like making them NaN (when the datatype is a floating point) is elegant and simplifies downstream code, but it means we have to be very careful not to call reduction functions (like min or sum) on them without masking out NaNs first.

[JohnAshburner]

I tend to use NaNs, and try to take care about handling them properly. Reduction functions are likely to give incorrect values whatever code is included in the data (and the 999 code for SPSS missing data has probably been the cause of a lot of dodgy stats). With NaNs, you at least see that you are not handling missing values properly. Julia has a framework for handling missingness. I don't know about Python/Torch.

[balbasty]

It should be done in this commit: 239060f9e90b76a7d525cb7ada0ac93b8fb8d421