Get a standard interface for the functions using the noise variance

[samuelstjean]

This is a bit confusing, but so far we have two method to estimate the noise, and almost no method to actually take advantage of what they offer nor point to the limitation of one or the other.

This should be made clearer, since people seems to get confused with their usage look at the mailing list for a recent discussion on restore and problems with dipy.denoise.estimate_sigma).

So, we should revamp the methods using incertainty in their formulation to accept global one valued input or local input for each voxel. So far, that would be restore and nlmeans, any other I'm not aware of/forgetting?

And also I suggest to make them more consistent. For example, restore has a built-in fix for rician noise computed from the variance of a gaussian distribution, but this does no make sense if you compute directly that value as piesno does. One way to uniformize this would be to put the fix inside the noise estimation module under a new function, and people could apply it themselves if need be.

Another thing is that no function accept a noise parameters on a voxelwise fashion, which is needed for spatially varying noise. It would be easier to have two codepath : one with a single value for the whole 4D volume, and one accepting a 3D volume for each voxel, thus covering all cases for all algorithm. And laos we should add in the doc that the noise_estimate function should really be used on a whole volume and not slices/subparts of the dataset, as it is based on a statistical estimator, i.e. you need a ton of data for it to be valid.

So that's it, I also have pretty picture to show why taking both single value and 3D volume of variance input is needed to cover all cases. So any other algorithm we should adapt while I'm at it?

[Garyfallidis]

Hi, this is an interesting proposition. I think nlmeans and restore are good candidates for now.

In your comment above you talk about having a picture showing the advantage of your suggestion but you are not showing us the picture.

[samuelstjean]

I was upon request of course ;) There it is, with a french caption. It basically says that the noise nature is ok if you stay at low spatial resolution, but becomes troublesome the moment you acquire higher spatial resolution data in diffusion MRI, since the signal is lower in each voxel and that the measure of interest is the signal loss.

So on the left a 1.2 mm acquistion, ans on the right same thing at 1.9 mm. For the left acquistion, you absolutely need to have a voxelwise estimation, while on the right picture you can probably manage fine with a single value for the whole volume. Both picture are of the same subject on the same scanner of course, only the voxelsize changed.

[arokem]

Yes - I think that would be an improvement.

For the time being, notice that the RESTORE implementation can accept an array of sigma values, if the user believes there is any reason that some directions might be noisier than others. This is not documented very well in the restore_fit_tensor function, but is explicit in the _nnls_err_func https://github.com/nipy/dipy/blob/master/dipy/reconst/dti.py#L1342. This input needs to be a 1D array with a single sigma value per volume, so the current functionality does not address the issue that you bring up.

Looking through the code, I actually see no particular block to extending the sigma kwarg to handle 3D volumes, in addition to the already existing support for 1D arrays (one sigma per volume) and for scalars (one global sigma). I think that you would need to introduce some logic into the TensorModel fit method to handle masking into this array, in a similar way that the signal is now handled, and indexing into this array in the restore_fit_tensor function. But I think that it's just local changes within these functions. With the obvious caveat that I have just looked it over, have not tried to actually implement it. Do you want to give it a go?

As for the images you showed: are these images of the noise or of the signal? If signal, is the non diffusion-weighted (b0) signal?

[samuelstjean]

For the easy part, these are both at b1000, bvec the closest to (0,1,0) if I remember correctly. You don't see that effect much on the b0 since the signal is really strong and relatively homogeneous, but here you go, another scan at 1mm, 1.25mm, 1.5mm and 1.9mm. You can actually see the signal drop with the voxel size (or b-value), while the noise is the same across the board.

For the second part of your comment, I'm not sure anything else than a single scalar or a 3D volume makes sense. What I mean by that is think as the noise as being sampled from a single magic distribution on a voxelwise basis. Disregarding the b-value or the sequence used, that distribution is sampled. Whenever the signal is high (T1, b0) the noise value is drowned in real signal, so no problem. Whenever you have a low signal (high b-value/high spatial resolution) you start to have more noise than signal, hence the impression that the noise is stronger, while it's more like the signal is getting weaker instead.

So a 1D array for a 4D volume does not make sense, as the noise caracteristic should be the same unless you have weird hardware problems (like spiking/dropout/whatever can go bad in complex space) I think.

As for the implementation details, the idea would be

1. Make nlmeans/RESTORE accept 1D/1 value or 3D arrays. We simply pick out the median from a whole 4D volume, which is stable for a single b-value, to get the sigma we use for nlmeans.
2. Specify that estimate_sigma is intended for whole volume estimation (remember it was developped for T1 noise estimation), and not neighborhood estimation since it breaks the statiscal property of the estimator as happening over at the mailing list I suspect.
3. Roll out a (real) 3D noise estimation procedure.

Regarding point 2, I suggest doc fixing to make it clearer. Point 1 is changing interface/patching code around, and point 3 is well, depends. I have two way already implemented for that, mine and the one from http://www.ncbi.nlm.nih.gov/pubmed/25499191 (they also give out matlab code). There is also one from http://www.ncbi.nlm.nih.gov/pubmed/25465845 which is available in R code and seems to work well, but takes so much time it's not realistically useable I think.

[arokem]

Thanks!

If I understand these images and your interpretation correctly, the noise
distribution is pretty much identical across voxels within each
acquisition, and across acquisitions, and it's just the signal amplitude
that actually changes. Is that correct? This seems consistent with our
findings comparing the noise in different b-values:
http://nbviewer.ipython.org/github/arokem/osmosis/blob/master/doc/paper_figures/AppendixB.ipynb
(TLDR; noise distrbution hardly changes for b=1,2,4k, while signal goes
down as approximately e^-b, no surprise). But, IIUC, this means that we can
use the same sigma for the entire volume. In the scenario where the noise
is the same in the entire volume (and only the signal changes), why would
we want to use a different sigma in each voxel?

I agree with you that using a different sigma for each volume is a strange
thing to do. This is a legacy that we inherited from @rfdougherty's Matlab
implementation, on which we based ours (if memory serves). If this confuses
people, we can refactor that out, or warn people that they are doing
something strange.

On Thu, Apr 30, 2015 at 10:54 AM, Samuel St-Jean notifications@github.com
wrote:

For the easy part, these are both at b1000, bvec the closest to (0,1,0) if
I remember correctly. You don't see that effect much on the b0 since the
signal is really strong and relatively homogeneous, but here you go,
another scan at 1mm, 1.25mm, 1.5mm and 1.9mm. You can actually see the
signal drop with the voxel size (or b-value), while the noise is the same
across the board.
[image: voxel_size]
https://cloud.githubusercontent.com/assets/3030760/7418705/d123908a-ef3d-11e4-869c-7e81576eb602.png

For the second part of your comment, I'm not sure anything else than a
single scalar or a 3D volume makes sense. What I mean by that is think as
the noise as being sampled from a single magic distribution on a voxelwise
basis. Disregarding the b-value or the sequence used, that distribution is
sampled. Whenever the signal is high (T1, b0) the noise value is drowned in
real signal, so no problem. Whenever you have a low signal (high
b-value/high spatial resolution) you start to have more noise than signal,
hence the impression that the noise is stronger, while it's more like the
signal is getting weaker instead.

So a 1D array for a 4D volume does not make sense, as the noise
caracteristic should be the same unless you have weird hardware problems
(like spiking/dropout/whatever can go bad in complex space) I think.

As for the implementation details, the idea would be

1. Make nlmeans/RESTORE accept 1D/1 value or 3D arrays. We simply pick
   out the median from a whole 4D volume, which is stable for a single
   b-value, to get the sigma we use for nlmeans.
2. Specify that estimate_sigma is intended for whole volume estimation
   (remember it was developped for T1 noise estimation), and not neighborhood
   estimation since it breaks the statiscal property of the estimator as
   happening over at the mailing list I suspect.
3. Roll out a (real) 3D noise estimation procedure.

Regarding point 2, I suggest doc fixing to make it clearer. Point 1 is
changing interface/patching code around, and point 3 is well, depends. I
have two way already implemented for that, mine and the one from
http://www.ncbi.nlm.nih.gov/pubmed/25499191 (they also give out matlab
code). There is also one from http://www.ncbi.nlm.nih.gov/pubmed/25465845
which is available in R code and seems to work well, but takes so much time
it's not realistically useable I think.

—
Reply to this email directly or view it on GitHub
#636 (comment).

[samuelstjean]

Yep, that it, noise should be the same, the signal changes though. And using partial fourier, acceleration (so always) changes the noise distribution, so the signal is lower on various spatial region AND the noise is also stronger, which amplifies the perceived effect. Now the next picture is actually pure noise (you simply dont measure the last echo in the diffusion sequence if I remember correctly) in a coronal view, and you can see its lower near the top since the antennas are there, while noisier near the neck.

[arokem]

I am convinced!

Let's make a plan:

From your list, number 2 is obviously the lowest-hanging fruit, because it
involves mostly the documentation. Do I understand correctly that you are
also suggesting we change the RESTORE example in the gallery? It seems that
we are using estimate_sigma incorrectly there?

Number 1 is not too hard, I think, so can be the next thing to do. Do you
want to give that a go?

Finally, number 3 as you say depends. When you say that you already have a
method implemented for that, are you referring to the PIESNO
implementation? Or is this something that you have implemented for your own
use, and is not ready to put in dipy yet? If the latter, it would be good
to take a look - maybe you could share this as a gist?

On Thu, Apr 30, 2015 at 11:28 AM, Samuel St-Jean notifications@github.com
wrote:

Yep, that it, noise should be the same, the signal changes though. And
using partial fourier, acceleration (so always) changes the noise
distribution, so the signal is lower on various spatial region AND the
noise is also stronger, which amplifies the perceived effect. Now the next
picture is actually pure noise (you simply dont measure the last echo in
the diffusion sequence if I remember ocrrectly) in a coronal view, and you
can see its lower near the top since the antennas are there, while noisier
near the neck.
[image: nmap]
https://cloud.githubusercontent.com/assets/3030760/7419738/276bb4b6-ef45-11e4-8a89-2039aeaab773.png

—
Reply to this email directly or view it on GitHub
#636 (comment).

[samuelstjean]

1. Is dead easy as it's pointing out flaws and warning doc wise.
2. Is also doable (and started) but for the restore example I don't know. Using one value should be the right thing to do I think (like taking the median on the 1D array, they should all be close for the b-values of the dataset), but it migth scare people away. It would be better to use 1 value and show the possible improvement with using a whole 3D array in the same example.
3. Is well, you have piesno which half does that since it assume the same variance inside a slice, but you get a 2D array since the value is different on each slice. The method I implemented is one of my own, it can use the noise map from the above image to estimate stuff or recreates one if you don't have it (because honestly, nobody except us is acquiring pure noise I think, and we only recently started doing that to test stuff around). And it's also unpublished (yet), based on nothing, etc. etc,

I have another implementation based on the first mentionned paper, but it only works on rician distribution by design, while my invention works on everything.

[samuelstjean]

Apparently I have a somewhat old but probably working version over here https://github.com/samuelstjean/dipy/blob/stabilisation/dipy/denoise/signal_transformation_framework.py

bottom completely, the last two functions.

[ShreyasFadnavis]

Thank you for this discussion! I guess we can close this issue for now. Feel free to reopen if anything on these lines pops up in the future. For diffusion data, we have MPPCA and Patch2Self both of which dont need an explicit sigma estimate!