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

Piesno #337

Merged
merged 18 commits into from May 30, 2014

Conversation

Projects
None yet
3 participants
@samuelstjean
Contributor

samuelstjean commented Mar 12, 2014

The is a followup of #277 in a new and improved version. This provides a way to estimate automagically the noise in 4D data, where the acquired data is the same amongst the last axis.

samuelstjean added some commits Mar 12, 2014

from os.path import join as pjoin, dirname
import nibabel as nib
test_piesno = nib.load(pjoin(dirname(__file__), 'test_piesno.nii.gz')).get_data()

This comment has been minimized.

@arokem

arokem Mar 12, 2014

Member

This is currently failing on travis (https://travis-ci.org/nipy/dipy/jobs/20637708#L3325). I wouldn't put this in the __init__.py file. Instead, put this in the test itself, if you really need it. Also, if you want to read files, use the schemes that exist in dipy.data. Is it possible to use one of the data-sets that is already defined there for this test? That would be preferable, I would think

This comment has been minimized.

@samuelstjean

samuelstjean Mar 12, 2014

Contributor

This one is given by the authors and they also give the value they get with their code, meaning we can check if this implementation gives exactly the same result. I could also run their version on any existing dataset, should both do the same thing.

This comment has been minimized.

@arokem

arokem Mar 12, 2014

Member

OK - one of two options, I think: either grab their data from their website
using a fetcher, or run their code on one of the data sets that we already
ship and benchmark to that.

On Wed, Mar 12, 2014 at 3:19 PM, Samuel St-Jean notifications@github.comwrote:

In dipy/denoise/tests/init.py:

@@ -0,0 +1,4 @@
+from os.path import join as pjoin, dirname
+import nibabel as nib
+
+test_piesno = nib.load(pjoin(dirname(file), 'test_piesno.nii.gz')).get_data()

This one is given by the authors and they also give the value they get
with their code, meaning we can check if this implementation gives exactly
the same result. I could also run their version on any existing dataset,
should both do the same thing.


Reply to this email directly or view it on GitHubhttps://github.com//pull/337/files#r10544335
.

This comment has been minimized.

@samuelstjean

samuelstjean Mar 12, 2014

Contributor

I'll go for the second idea, since the data is jammed in a big zip file
full of stuff.
Le 2014-03-12 18:21, Ariel Rokem a écrit :

In dipy/denoise/tests/init.py:

@@ -0,0 +1,4 @@
+from os.path import join as pjoin, dirname
+import nibabel as nib
+
+test_piesno = nib.load(pjoin(dirname(file), 'test_piesno.nii.gz')).get_data()
OK - one of two options, I think: either grab their data from their
website using a fetcher, or run their code on one of the data sets
that we already ship and benchmark to that. On Wed, Mar 12, 2014 at
3:19 PM, Samuel St-Jean notifications@github.comwrote:
… <#>
In dipy/denoise/tests/init.py: > @@ -0,0 +1,4 @@ > +from os.path
import join as pjoin, dirname > +import nibabel as nib > + >
+test_piesno = nib.load(pjoin(dirname(file),
'test_piesno.nii.gz')).get_data() This one is given by the authors and
they also give the value they get with their code, meaning we can
check if this implementation gives exactly the same result. I could
also run their version on any existing dataset, should both do the
same thing. — Reply to this email directly or view it on
GitHubhttps://github.com//pull/337/files#r10544335 .


Reply to this email directly or view it on GitHub
https://github.com/nipy/dipy/pull/337/files#r10544432.

This comment has been minimized.

@arokem

arokem Mar 12, 2014

Member

Awesome

On Wed, Mar 12, 2014 at 3:23 PM, Samuel St-Jean notifications@github.comwrote:

In dipy/denoise/tests/init.py:

@@ -0,0 +1,4 @@
+from os.path import join as pjoin, dirname
+import nibabel as nib
+
+test_piesno = nib.load(pjoin(dirname(file), 'test_piesno.nii.gz')).get_data()

I'll go for the second idea, since the data is jammed in a big zip file
full of stuff. Le 2014-03-12 18:21, Ariel Rokem a écrit :
… <#144b8647ef430658_>
In dipy/denoise/tests/init.py: > @@ -0,0 +1,4 @@ > +from os.path
import join as pjoin, dirname > +import nibabel as nib > + > +test_piesno =
nib.load(pjoin(dirname(file), 'test_piesno.nii.gz')).get_data() OK -
one of two options, I think: either grab their data from their website
using a fetcher, or run their code on one of the data sets that we already
ship and benchmark to that. On Wed, Mar 12, 2014 at 3:19 PM, Samuel St-Jean
notifications@github.comwrote: … <#> In dipy/denoise/tests/init.py:

@@ -0,0 +1,4 @@ > +from os.path import join as pjoin, dirname > +import
nibabel as nib > + > +test_piesno = nib.load(pjoin(dirname(file),
'test_piesno.nii.gz')).get_data() This one is given by the authors and they
also give the value they get with their code, meaning we can check if this
implementation gives exactly the same result. I could also run their
version on any existing dataset, should both do the same thing. — Reply to
this email directly or view it on GitHub<
https://github.com/nipy/dipy/pull/337/files#r10544335> . — Reply to this
email directly or view it on GitHub <
https://github.com/nipy/dipy/pull/337/files#r10544432>.


Reply to this email directly or view it on GitHubhttps://github.com//pull/337/files#r10544520
.

Parameters
-----------
data : numpy array

This comment has been minimized.

@Garyfallidis

Garyfallidis Mar 17, 2014

Member

write here: array or ndarray

data : numpy array
The magnitude signals to analyse. The last dimension must contain the
same realisation of the volume, such as dMRI or fMRI data.

This comment has been minimized.

@Garyfallidis

Garyfallidis Mar 17, 2014

Member

What do you you mean here?

This comment has been minimized.

@samuelstjean

samuelstjean Mar 17, 2014

Contributor

That it will work on repeated measurement data, but not on structural data. That means basically dMRI of fMRI, but not T1/T2/etc.

This comment has been minimized.

@Garyfallidis
def piesno(data, N=1, alpha=0.01, l=100, itermax=100, eps=1e-10):
"""
A routine for finding the underlying gaussian distribution standard
deviation from magnitude signals.

This comment has been minimized.

@Garyfallidis

Garyfallidis Mar 17, 2014

Member

What do you mean magnitude signals? Aren't these signals what we always use in dMRI?

This comment has been minimized.

@samuelstjean

samuelstjean Mar 17, 2014

Contributor

Usually, but you have a slew of other reconstruction methods that might not give that, but they are probably seldom used.

The magnitude signals to analyse. The last dimension must contain the
same realisation of the volume, such as dMRI or fMRI data.
N : int

This comment has been minimized.

@Garyfallidis

Garyfallidis Mar 17, 2014

Member

How important is these parameter? Often people may not know the number of phase coils.

This comment has been minimized.

@samuelstjean

samuelstjean Mar 17, 2014

Contributor

It defines the probability distribution of the whole acquisition, so it's needed to properly find the real value of sigma. It's also an hassle to find (we don't even know what is supposed to be the real value at the CHUS), so people will probably just put N=1 for rician noise.

This comment has been minimized.

@samuelstjean

samuelstjean Mar 17, 2014

Contributor

Does anyone know if this information is supposed to be in the header or something? Couldn't find anything about that either.

def piesno(data, N=1, alpha=0.01, l=100, itermax=100, eps=1e-10):
"""
A routine for finding the underlying gaussian distribution standard

This comment has been minimized.

@Garyfallidis

Garyfallidis Mar 17, 2014

Member

Can you say here what PIESNO stands for?

import numpy as np
from numpy.testing import assert_almost_equal
from dipy.denoise.signal_transformation_framework import _inv_nchi_cdf, piesno

This comment has been minimized.

@Garyfallidis

Garyfallidis Mar 17, 2014

Member

This looks wrong... there is no signal_transformation_framework file. Right?

# Values taken from hispeed.MedianPIESNO with the test data
# in the package computed in matlab
sigma = piesno(test_piesno, N=8, alpha=0.01, l=1)

This comment has been minimized.

@Garyfallidis

Garyfallidis Mar 17, 2014

Member

Okay, your testing needs reworking here. You use test_piesno as a parameter here but that is nowhere to find.

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Mar 17, 2014

Please also resolve the problem with the datasets that you use for the tests. For tests we avoid fetching online data and try to use either small data from dipy.data or create synthetic data on the fly. However, in the tutorial you should have some real data that you can fetch.

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Mar 17, 2014

Information about the coils is kept in the dicom headers but it is not available in the niftis.

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Apr 15, 2014

Hi Sam, sorry for the delay but can you rebase this so we can see the coverage of your tests with Travis?

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Apr 15, 2014

Should be good I think. Discovered that this could be used by the RESTORE algorithm to estimate sigma robustly btw.

As for code coverage, there is also https://coveralls.io/ that seems popular.

@arokem

This comment has been minimized.

Member

arokem commented Apr 15, 2014

If you feel like updating the RESTORE example with this, that would be
excellent.

On Mon, Apr 14, 2014 at 5:36 PM, Samuel St-Jean notifications@github.comwrote:

Should be good I think. Discovered that this could be used by the RESTORE
algorithm to estimate sigma robustly btw.


Reply to this email directly or view it on GitHubhttps://github.com//pull/337#issuecomment-40434029
.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Apr 15, 2014

I'll probably update the SNR estimation, nlmeans and RESTORE examples at
the same time, in another PR. It would be cleaner this way.

Le 2014-04-15 11:08, Ariel Rokem a écrit :

If you feel like updating the RESTORE example with this, that would be
excellent.

On Mon, Apr 14, 2014 at 5:36 PM, Samuel St-Jean
notifications@github.comwrote:

Should be good I think. Discovered that this could be used by the
RESTORE
algorithm to estimate sigma robustly btw.


Reply to this email directly or view it on
GitHubhttps://github.com//pull/337#issuecomment-40434029
.


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

# Initial estimation of sigma
denom = np.sqrt(2 * _inv_nchi_cdf(N, 1, q))
m = np.percentile(data, q*100) / denom

This comment has been minimized.

@Garyfallidis
m = np.percentile(data, q*100) / denom
phi = np.arange(1, l + 1) * m/l
K = data.shape[-1]
sum_m2 = np.sum(data**2, axis=-1)

This comment has been minimized.

@Garyfallidis

This comment has been minimized.

@samuelstjean

samuelstjean May 27, 2014

Contributor

Having the exponent float with a space seems ugly to me.

This comment has been minimized.

@arokem

arokem May 27, 2014

Member

You mean complying to PEP8 in this line would be ugly? I agree with you. I
think we're usually reasonable about these things. In other words, no need
to put the ugly spaces in here. Do you have any sense for why the test is
currently failing on Travis?

https://travis-ci.org/nipy/dipy/jobs/26140067

On Tue, May 27, 2014 at 8:25 AM, Samuel St-Jean notifications@github.comwrote:

In dipy/denoise/noise_estimate.py:

  •               16: 0.5900487123737876,
    
  •               32: 0.5641772300866416,
    
  •               64: 0.5455611840489607,
    
  •              128: 0.5322811923303339}
    
  • if N in opt_quantile:
  •    q = opt_quantile[N]
    
  • else:
  •    q = 0.5
    
  • Initial estimation of sigma

  • denom = np.sqrt(2 * _inv_nchi_cdf(N, 1, q))
  • m = np.percentile(data, q*100) / denom
  • phi = np.arange(1, l + 1) * m/l
  • K = data.shape[-1]
  • sum_m2 = np.sum(data**2, axis=-1)

Having the exponent float with a space seems ugly to me.


Reply to this email directly or view it on GitHubhttps://github.com//pull/337/files#r13083485
.

if np.abs(sig - sig_prev) < eps:
break
s = sum_m2 / (2*K*sig**2)

This comment has been minimized.

@Garyfallidis
sig_prev = sig
# Numpy percentile must range in 0 to 100, hence q*100
sig = np.percentile(omega, q*100) / denom

This comment has been minimized.

@Garyfallidis
def test_inv_nchi():
# Values taken from hispeed.MedianPIESNO.lambdaPlus

This comment has been minimized.

@Garyfallidis

Garyfallidis Apr 18, 2014

Member

What are these?

This comment has been minimized.

@samuelstjean

samuelstjean Apr 18, 2014

Contributor

The reference, closed source java implementation from which I took the values.

This comment has been minimized.

@Garyfallidis

Garyfallidis Apr 18, 2014

Member

Okay add a link please in the comments.

This comment has been minimized.

@samuelstjean

samuelstjean Apr 18, 2014

Contributor

we'll, you actually need ot write an email to the guy to get it and sign a small disclaimer, and the values are in the paper. Much simpler for anyone to just check that against the paper version I think.

# in the package computed in matlab
test_piesno_data = nib.load(dipy.data.get_data("test_piesno")).get_data()
sigma = piesno(test_piesno_data, N=8, alpha=0.01, l=1)[0]
assert_almost_equal(sigma, 0.010749458025559)

This comment has been minimized.

@Garyfallidis

Garyfallidis Apr 18, 2014

Member

What are you actually trying to test here? Why the value should be this one?

This comment has been minimized.

@samuelstjean

samuelstjean Apr 18, 2014

Contributor

They are given in the paper for the dataset that comes with the examples.

This comment has been minimized.

@Garyfallidis

Garyfallidis Apr 18, 2014

Member

Okay write in the comment the page of the paper that you found these values.

The estimated standard deviation of the gaussian noise.
mask : ndarray
A boolean mask indicating the voxels identified as pure noise.

This comment has been minimized.

@Garyfallidis

Garyfallidis Apr 18, 2014

Member

Can you make the mask optional?

References
------------
.. [1]. Koay CG, Ozarslan E and Pierpaoli C.

This comment has been minimized.

@Garyfallidis

Garyfallidis May 28, 2014

Member

Typo: [1]. should be [1]

Parameters
-----------

This comment has been minimized.

@Garyfallidis

Garyfallidis May 28, 2014

Member

Remove space between Parameters and data

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented May 30, 2014

Thx! We now need to update the denoising and restore tutorials to use PIESNO.

Garyfallidis added a commit that referenced this pull request May 30, 2014

@Garyfallidis Garyfallidis merged commit aa65931 into nipy:master May 30, 2014

1 check passed

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

@samuelstjean samuelstjean deleted the samuelstjean:piesno branch Jul 4, 2014

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