New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added local variance estimation #433

Merged
merged 16 commits into from Nov 13, 2014

Conversation

Projects
None yet
5 participants
@samuelstjean
Contributor

samuelstjean commented Oct 3, 2014

Fixes #432, Seems to gives the same thing as our closed source implementation of nlmeans so far.

This PR :

In [109]: estimate_sigma(d)
Out[109]:
array([ 31.94007111, 13.55418587, 13.74776077, 13.63606453,
13.67319965, 13.60422993, 13.56204987, 13.60548496,
13.76000309, 13.57973385, 13.72612858, 13.7350111 ,
13.65256119, 13.58515263, 13.66841412, 13.60126591,
13.58483315, 13.58515263, 13.643466 , 13.74186325,
13.59409237, 13.73767948, 13.60761738, 13.59641361,
13.73959732, 13.59790134, 13.73537064, 13.66113281,
13.63371849, 13.6424818 , 13.60521698, 13.59016132,
13.59904766, 13.62820435, 13.5684433 , 13.56952667,
13.57829094, 13.75860214, 13.63626385, 13.6825037 ,
13.62080669, 13.65673733, 13.69610214, 13.63846493,
13.74402618, 13.71609783, 13.57646942, 13.70084381,
13.6728611 , 13.65078163, 13.6009016 , 13.70344543,
13.71196938, 13.72620392, 13.63866806, 13.67575741,
13.74444294, 13.64913177, 13.7595911 , 13.60122585,
13.76624775, 13.66040325, 13.64105892, 13.72773266, 13.8091774 ], dtype=float32)

The other implementation :
Estimated standard deviation of noise :31.9401 --> Sigma = 31.9401
Estimated standard deviation of noise :13.5542 --> Sigma = 13.5542

I alos added a cheap check to remove zeros when the background is absent, which won't give the same numbers, but is close on the same test (the data has some background set to zero)

In [94]: estimate_sigma(d)
Out[94]:
array([ 32.0320015 , 13.59815121, 13.79326344, 13.68167877,
13.71766853, 13.64865971, 13.60599041, 13.64930344,
13.80526352, 13.62382603, 13.77113152, 13.77896404,
13.69764805, 13.62915134, 13.71269417, 13.64586449,
13.6284008 , 13.62946415, 13.68790531, 13.78582764,
13.63778496, 13.78283215, 13.65133572, 13.63971615,
13.78512955, 13.64159966, 13.77961826, 13.70498466,
13.67730331, 13.68666077, 13.64945412, 13.63408661,
13.64239788, 13.67206287, 13.61290073, 13.61320686,
13.62299824, 13.80318546, 13.68140888, 13.72596836,
13.66493225, 13.70094013, 13.74094391, 13.68257999,
13.78874302, 13.75980568, 13.62131596, 13.74606228,
13.71796322, 13.69530106, 13.64515877, 13.74870014,
13.75715065, 13.77051258, 13.68297386, 13.7193985 ,
13.78968048, 13.69384766, 13.80404091, 13.64589691,
13.81126785, 13.70522404, 13.68564224, 13.77248096, 13.85372639], dtype=float32)

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Oct 3, 2014

Good. Tests please!

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Oct 3, 2014

Ok, test what... Best we can do is compute something and check I did the same bug as the reference implementation, and it only gives 4 significant numbers.

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Oct 3, 2014

Create a simulated image where you can calculate the exact value of the sigma. And then test that you get that value.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Oct 3, 2014

Added fancy tests.

samuelstjean added some commits Oct 7, 2014

@@ -0,0 +1,67 @@
from __future__ import division, print_function

This comment has been minimized.

@Garyfallidis

Garyfallidis Oct 11, 2014

Member

There must be something wrong with the name of this file.

This comment has been minimized.

@Garyfallidis

Garyfallidis Oct 11, 2014

Member

You must mean test_noise_estimate.py . Correct please...

samuelstjean added some commits Oct 11, 2014

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Nov 2, 2014

Hi @samuelstjean,

There is clearly a problem. In the following code I am using the T1 image provided with the Stanford data. So, you should have it too. Check your .dipy folder or run the connectivity tutorial to fetch it.

import nibabel as nib
from dipy.denoise.noise_estimate import estimate_sigma
from dipy.denoise.nlmeans import nlmeans

fname = '/home/eleftherios/Desktop/t1.nii.gz'
fname2 = '/home/eleftherios/Desktop/t1_denoised.nii.gz'

img = nib.load(fname)
data = img.get_data()
affine = img.get_affine()

sigma = estimate_sigma(data)
data2 = nlmeans(data, sigma[0])

nib.save(nib.Nifti1Image(data2, affine), fname2)

The resulted denoised image is extremely over smoothed. Try it for yourself.
What is the problem?

@arokem

This comment has been minimized.

Member

arokem commented Nov 2, 2014

Hey @Garyfallidis - just a note - you probably got the T1 when you reviewed
the LiFE PR. I added that as a data-set there. @samulelstjean : you can get
the T1 file here: http://purl.stanford.edu/yx282xq2090/

On Sun, Nov 2, 2014 at 3:55 PM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

Hi @samuelstjean https://github.com/samuelstjean,

There is clearly a problem. In the following code I am using the T1 image
provided with the Stanford data. So, you should have it too. Check your
.dipy folder or run the connectivity tutorial to fetch it.

import nibabel as nib
from dipy.denoise.noise_estimate import estimate_sigma
from dipy.denoise.nlmeans import nlmeans

fname = '/home/eleftherios/Desktop/t1.nii.gz'
fname2 = '/home/eleftherios/Desktop/t1_denoised.nii.gz'

img = nib.load(fname)
data = img.get_data()
affine = img.get_affine()

sigma = estimate_sigma(data)
data2 = nlmeans(data, sigma[0])

nib.save(nib.Nifti1Image(data2, affine), fname2)

The resulted denoised image is extremely over smoothed. Try it for
yourself.
What is the problem?


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

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 3, 2014

The problem is in the dataset itself, or rather the noise estimation performs poorly on it. Seems to be a datatype issue since I get a different answer

In [16]: estimate_sigma(data.astype(np.float32), False)
Out[16]: array([ 3755.09228516], dtype=float32)

In [17]: estimate_sigma(data.astype(np.int16), False)
Out[17]: array([ 13012.04296875], dtype=float32)

Of course the float32 answer is more in agreement with our old implementation, which gives
Estimated standard deviation of noise :3755.09 --> Sigma = 3755.09

So the culprit could be the convolve from scipy, seems like it's the only potential problematic line. The function doesn't support specifying a return dtype, but we can't convert everything to float for memory usage sake either. Any suggestion?

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 3, 2014

And as a side note, the T1 is not super high quality in my opinion. Yes it goes hand in hand with the stanford dataset, but the spatial resolution is poor to be used in a standalone example at 2mm^3.

@arokem

This comment has been minimized.

Member

arokem commented Nov 3, 2014

Agreed. I can make the original data (at 0.9 mm isotropic) also available.
The tricky bit is that to resample and align from one to the other, I used
functions from nipy, and I think we would want to not have that be a
dependency, so I just uploaded the result of the resampling for now. What
do you think?

On Sun, Nov 2, 2014 at 4:20 PM, Samuel St-Jean notifications@github.com
wrote:

And as a side note, the T1 is not super high quality in my opinion. Yes it
goes hand in hand with the stanford dataset, but the spatial resolution is
poor to be used in a standalone example at 2mm^3.


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

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 3, 2014

Well, you could make both available, I guess you resampled the T1 to the same size as the DWIs? That way, one could either use the downsampled T1, or upsample the diffusion to the T1 space (which one would need to do himself since it goes out of the scope of the example).

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Nov 3, 2014

I don't think that you should resample the data to help the noise estimation. The noise estimation should work with this dataset and be uneffected of dtype or similar problems.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 3, 2014

I never said that, I said if you want to use it in an example, another T1 might be better.

As for the datatype thing, I found the problem, it's really the dataset actually.

In [22]: estimate_sigma(data.astype(np.int32), False)
Out[22]: array([ 3755.09204102], dtype=float32)

In [16]: data.min()
Out[16]: 0

In [17]: data.max()
Out[17]: 32767

it's at the limit of 2**15, so it's looping around probably inside the convolve.

@rfdougherty

This comment has been minimized.

rfdougherty commented Nov 3, 2014

The convolve is probably overflowing the voxels at or near 32767, thus the inflated sigma estimates. This seems like a limitation in convolve-- it doesn't clamp results to dtype max/min and thus allows overflow. It's fortunate that this dataset had some extreme values and thus revealed the issue. 32767 is a perfectly valid value for an int16 image, so the code should not fail with such valid input.

Converting to a float where necessary seems reasonable to me. If done in the loop, it won't increase memory usage much for most real-world cases. You might be able to do it by specifying an output for convolve rather than converting the array to float. E.g.:

conv_out = np.zeros(arr[...,0].shape, dtype='float64')
...
convolve(arr[...,i], k, output=conv_out)
mean_block = np.sqrt(6/7) * (arr[..., i] - 1/6 * conv_out)

The only other option I can see is to modify the input data, clipping it to dtype max - 6 and dtype min + 6 (unless dtype min==0).

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 6, 2014

Now that we have a cheap global estimation technique, should we try a cheapo local one also to cope with non stationnary noise? This could turn the NLMEANS implementation into an adaptive NLMEANS instead.

@mdesco

This comment has been minimized.

Contributor

mdesco commented Nov 6, 2014

On Wed, Nov 5, 2014 at 8:52 PM, Samuel St-Jean notifications@github.com
wrote:

Now that we have a cheap global estimation technique, should we try a
cheapo local one also to cope with non stationnary noise?

Yes we should.

This could turn the NLMEANS implementation into an adaptive NLMEANS
instead.

That would be a nice improvement.


Reply to this email directly or view it on GitHub
https://github.com/nipy/dipy/pull/433#issuecomment-61915531.[image: Web
Bug from
https://github.com/notifications/beacon/ACEgB1Haihm80l1LvMVlkzGlFG0NOEKvks5nKsvegaJpZM4CqpN6.gif]

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Nov 13, 2014

Hi @samuelstjean,
It seems it is working now :). Thank you @rfdougherty for the feedback.

But what is strange with function estimate_sigma is that you return a list even if you have one volume only i.e. one sigma. On the other hand in nlmeans you don't support multiple sigmas for different volumes. Can you elaborate on that?

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 13, 2014

It does, look at line 40 and 43 in nlmeans.py. The 4D version iterates on a list of sigma values now

sigma_arr = np.ones(arr.shape[-1], dtype=np.float32) * sigma

If you give an array for each volume or a single value it will work.

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Nov 13, 2014

Okay and what about returning a float rather than a list when there is only one volume? Can you correct that?

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 13, 2014

I fear that would break the generality of the code, and a single item list can be used in the same way as a float right?

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Nov 13, 2014

Okay this was just a suggestion. Your call!

Garyfallidis added a commit that referenced this pull request Nov 13, 2014

Merge pull request #433 from samuelstjean/local_variance
Added local variance estimation

@Garyfallidis Garyfallidis merged commit ac24533 into nipy:master Nov 13, 2014

1 check passed

continuous-integration/travis-ci The Travis CI build passed
Details

@samuelstjean samuelstjean deleted the samuelstjean:local_variance branch Nov 13, 2014

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Nov 13, 2014

Merged! Congratulations Sam!
@arokem can you update asap the RESTORE tutorial and @samuelstjean can you update the nlmeans tutorial?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment