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 example #576

Merged
merged 23 commits into from Mar 29, 2015

Conversation

Projects
None yet
5 participants
@arokem
Member

arokem commented Feb 19, 2015

This supersedes #390

for a in ax:
a.set_axis_off()
# Uncomment the coming line if you want a window to show the images.

This comment has been minimized.

@arokem

arokem Feb 19, 2015

Member

Let's take this out...

"""
.. [Koay2009] Koay C.G., E. Ozarslan, C. Pierpaoli. Probabilistic Identification and Estimation of Noise (PIESNO): A self-consistent approach and its applications in MRI. JMR, 199(1):94-103, 2009.

This comment has been minimized.

@arokem

arokem Feb 19, 2015

Member

This line needs to be wrapped

sigma[num] = sig
mask[num] = idx
if return_mask:

This comment has been minimized.

@samuelstjean

samuelstjean Feb 19, 2015

Contributor

does this PR actually remove the old noise estimation stuff? It's still useful for T1 imaging if that's the case.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Feb 19, 2015

If there is a way to get those, my fixes for all the previous comment are actually over there mdesco#10

@arokem

This comment has been minimized.

Member

arokem commented Feb 19, 2015

Best would be to rebase your branch on master and then make a new PR
against my recent PR (#576). I promise to
be quick about merging in your proposed changes, so we can keep going.

On Thu, Feb 19, 2015 at 12:59 PM, Samuel St-Jean notifications@github.com
wrote:

If there is a way to get those, my fixes for all the previous comment are
actually over there mdesco#10 mdesco#10


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

@arokem arokem force-pushed the arokem:mdesco-piesno-example branch from 29f7ce5 to 42c1308 Feb 22, 2015

@arokem arokem referenced this pull request Feb 22, 2015

Closed

Piesno fix #12

@arokem arokem force-pushed the arokem:mdesco-piesno-example branch from 5426f02 to b52e5cb Feb 23, 2015

@arokem arokem referenced this pull request Feb 23, 2015

Merged

Mdesco piesno example #13

import numpy as np
from scipy.special import gammainccinv
from scipy.ndimage.filters import convolve
from scipy.stats import mode

This comment has been minimized.

@samuelstjean

samuelstjean Feb 23, 2015

Contributor

the mode import is not used anymore, so maybe some other random stuff got in again :/

This comment has been minimized.

@arokem

arokem Feb 23, 2015

Member

No - that's my bad. Wasn't sure what to do with this line in the rebase.
Keep looking, though. I might've messed up in other places...

On Mon, Feb 23, 2015 at 12:48 PM, Samuel St-Jean notifications@github.com
wrote:

In dipy/denoise/noise_estimate.py
#576 (comment):

import numpy as np

from scipy.special import gammainccinv
from scipy.ndimage.filters import convolve
+from scipy.stats import mode

the mode import is not used anymore, so maybe some other random stuff got
in again :/


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

# if the percentile is 0, then more than half of the slice is zero and give
# up:
if m == 0:

This comment has been minimized.

@samuelstjean

samuelstjean Feb 24, 2015

Contributor

Only thing missing if if we really want to have this check, or instead have a sanity check for when the whole slice is zero (and not only the quantile).

This comment has been minimized.

@arokem

arokem Feb 24, 2015

Member

Sorry, I am not sure I understand.

These lines are identical in this PR, and in the PR you did against this
branch, here:

arokem#13

Do you have a further suggestion for change? Should we change this line to
the following?

if m==0 and np.all(data==0)

On Tue, Feb 24, 2015 at 1:47 PM, Samuel St-Jean notifications@github.com
wrote:

In dipy/denoise/noise_estimate.py
#576 (comment):

 # Initial estimation of sigma
 denom = np.sqrt(2 * _inv_nchi_cdf(N, 1, q))
 m = np.percentile(data, q * 100) / denom
  • if the percentile is 0, then more than half of the slice is zero and give

  • up:

  • if m == 0:

Only thing missing if if we really want to have this check, or instead
have a sanity check for when the whole slice is zero (and not only the
quantile).


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

This comment has been minimized.

@samuelstjean

samuelstjean Feb 24, 2015

Contributor

It' an open suggestion let's say. I was playing with a very noisy, background masked dataset, and it cut 2-3 lines of gm/wm at the top of the coronal slice for this reason I suspect. Relaxing this warning for np.all(data==0) would probably work, since the algorithm already exits if no noise voxel is found , but it's computed if we decide to remove this check + warning.

This comment has been minimized.

@arokem

arokem Feb 26, 2015

Member

So you are actually suggesting to replace it by the following?

np.all(data==0)

This comment has been minimized.

@samuelstjean

samuelstjean Feb 26, 2015

Contributor

yep, that would be more mathematically sound to me. I don't really remember why people suggested that actually, maybe they found a weird corner case on real data or something, but so far I always get a plausible mask except on cropped datasets where the first quantile estimation is screwed up to be high instead. Maybe that's actually what they meant... So far, as long as there is background (even zeros), it works well. Crop the dataset, and you screw the maths also.

This comment has been minimized.

@arokem

arokem Feb 26, 2015

Member

Maybe we can also get someone else's opinion on this? @mdesco - does this make sense?

Replacing if m==0: with if np.all(data==0)?

This comment has been minimized.

@samuelstjean

samuelstjean Feb 26, 2015

Contributor

We are at the MRI for the full day tomorrow, I can try asking him in between scans.

This comment has been minimized.

@arokem

arokem Feb 26, 2015

Member

Cool. Please do. Happy scanning!

On Wed, Feb 25, 2015 at 6:17 PM, Samuel St-Jean notifications@github.com
wrote:

In dipy/denoise/noise_estimate.py
#576 (comment):

 # Initial estimation of sigma
 denom = np.sqrt(2 * _inv_nchi_cdf(N, 1, q))
 m = np.percentile(data, q * 100) / denom
  • if the percentile is 0, then more than half of the slice is zero and give

  • up:

  • if m == 0:

We are at the MRI for the full day tomorrow, I can try asking him in
between scans.


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

This comment has been minimized.

@samuelstjean

samuelstjean Feb 27, 2015

Contributor

completely forgot about it, our day was crammed already, , oh well he will probably see the message nevertheless

@arokem arokem force-pushed the arokem:mdesco-piesno-example branch from bf9c42b to 8379f6f Mar 2, 2015

@stefanv stefanv modified the milestone: 0.9 Mar 4, 2015

def _piesno_3D(data, N, alpha=0.01, l=100, itermax=100, eps=1e-5,
return_mask=False):

This comment has been minimized.

@fmorency

fmorency Mar 19, 2015

The return_mask parameter is always True (see above).

This comment has been minimized.

@samuelstjean

samuelstjean Mar 19, 2015

Contributor

Yes, the wrapper always calls the inner function with return_mask=True, then discards the mask if the user did not want it in the end at line 110. The 3D wrapped function does not return it, since it confused some people who mistook that as a brain mask to pass to a reconstruction algorithm for example, while it is instead a noise mask, identifying noisy voxels.

This comment has been minimized.

@arokem

arokem Mar 19, 2015

Member

Yes - I think this is fine. This is a 'private' function anyhow, so users are discouraged from using it directly.

reached if two subsequent estimates are smaller than eps.
Default: 1e-5.
return_mask : bool (optiona)

This comment has been minimized.

@fmorency

fmorency Mar 19, 2015

Typo. *optional

of 3 coil elements each. The signal is therefore received through 4 distinct
groups of receiver channels, yielding N = 4. Had we used a GE acquisition, we
would have used N=1 even if multiple channel coils are used because GE uses a
SENSE acceleration, which has a Rician noise nature and thus N is always 1.

This comment has been minimized.

@fmorency

fmorency Mar 19, 2015

SENSE reconstruction

As a convenience, we will estimate the noise for the whole volume in one go,
but it is also possible to get a slice by slice estimation of the noise if
it is more desirable through the :func:`dipy.denoise.noise_estimate.piesno_3D`

This comment has been minimized.

@fmorency

fmorency Mar 19, 2015

The method dipy.denoise.noise_estimate.piesno_3D doesn't exist anymore.

This comment has been minimized.

@samuelstjean

samuelstjean Mar 19, 2015

Contributor

Good catch, that example need some reworking since the function would return one single value at some point to not confuse people since no function at the time actually supported that, which was later changed to work as intended, and instead patch nlmeans to now accept arrays as input.

Anyhow, one can still get a single value from the function if required (by taking the mode or the median for example), but it's not hidden internally anymore.

This comment has been minimized.

@arokem

arokem Mar 19, 2015

Member

OK - I removed this paragraph. It is really obscure right now.

@arokem arokem force-pushed the arokem:mdesco-piesno-example branch from 8379f6f to 2e27209 Mar 19, 2015

@arokem

This comment has been minimized.

Member

arokem commented Mar 19, 2015

Thanks for the friendly nudge, @fmorency, and thanks for the review. If you think this is now ready to go, I will merge this into master (in a couple of days, so that others can chime in here), so that you can easily use this. I will ask you (and others) to keep an eye out, and report back with any problems/issues that you encounter.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Mar 19, 2015

Before that I would remove the check around line 200 to look for explicitly the sane value, as opposed to a threshold of zeros.

@arokem

This comment has been minimized.

Member

arokem commented Mar 19, 2015

Like so? arokem@00d7007

@arokem

This comment has been minimized.

Member

arokem commented Mar 19, 2015

No - that would probably not do. Let me rebase that out. What is the "sane value" here and how would you find it? Specifics, please.

@arokem arokem force-pushed the arokem:mdesco-piesno-example branch from 00d7007 to 2956ec9 Mar 19, 2015

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Mar 19, 2015

That would be having an entire slice made of the same value. Remember there is no sane default, we are hacking togegther a convient wrapper for people, which is not covered in the original paper at all. The only in doing this is I don't really know anymore. The saved computation time is probably meaningless, since the function discards invalid estimate and exits when there are none, so it works fine without such a hack.

@arokem

This comment has been minimized.

Member

arokem commented Mar 19, 2015

So, replace the current check for m==0 with the following check?

if np.var(slice) < tol:
  ...

On Thu, Mar 19, 2015 at 12:17 PM, Samuel St-Jean notifications@github.com
wrote:

That would be having an entire slice made of the same value. Remember
there is no sane default, we are hacking togegther a convient wrapper for
people, which is not covered in the original paper at all. The only in
doing this is I don't really know anymore. The saved computation time is
probably meaningless, since the function discards invalid estimate and
exits when there are none, so it works fine without such a hack.


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

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Mar 19, 2015

actually, juste removing that check altogether would probably still work fine. I've yet to find a real dataset to throw at it where it fails.

The only time I got the function to blatantly lie to me is on synthetic data when there is more data than background, then the first estimate of the quantile is screwed and the function reversely picks the data instead of the background. Of course this is super obvious while looking at the mask since everything is inverted.

The only way to get that on a real dataset would be to mask it and then crop it to the minimal bounding box, in which case the hypothesis for the function is void and I have no idea how to find out this case.

That would be the ideal intention behind those lines in all honesty, which is not really covered currently.

@arokem

This comment has been minimized.

Member

arokem commented Mar 19, 2015

OK - it's actually also explicitly stated in the Notes section of the
documentation, and should be obvious to users if they read the paper and
the example.

How's this?

arokem@bd11490

On Thu, Mar 19, 2015 at 12:33 PM, Samuel St-Jean notifications@github.com
wrote:

actually, juste removing that check altogether would probably still work
fine. I've yet to find a real dataset to throw at it where it fails.

The only time I got the function to blatantly lie to me is on synthetic
data when there is more data than background, then the first estimate of
the quantile is screwed and the function reversely picks the data instead
of the background. Of course this is super obvious while looking at the
mask since everything is inverted.

The only way to get that on a real dataset would be to mask it and then
crop it to the minimal bounding box, in which case the hypothesis for the
function is void and I have no idea how to find out this case.

That would be the ideal intention behind those lines in all honesty, which
is not really covered currently.


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

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Mar 19, 2015

Well, we are back to the original version, which is fine by me actually, but some people said that we should not assume that people read the doc, so I'm at a loss at how to cover this case.

Anyway, it's still better with all these usability notes since the current version didn't not explicitely mention one needs to wrap this function, leading to troubles in using it for everyone. So wither we assume people will read those notes, or we try to hack a way to find the case where the data is cropped.

@arokem

This comment has been minimized.

Member

arokem commented Mar 27, 2015

I was just waiting for you to test this out, and get your OK to merge this.

I am not sure what you want to address with additional work on this PR. If
you want to add anything, maybe do that on a separate PR, after we merge
this one? As mentioned earlier in this thread, the version of PIESNO on
current master is broken, and this fixes it.

On Fri, Mar 27, 2015 at 8:05 AM, Samuel St-Jean notifications@github.com
wrote:

Ok, since nobody seems to want to add anything, I'll think I'll go for a
wrapper + private function that always return the mask. github is dying
right now, so I'll make a PR agaisn't this branch with those changes and we
can see what people think about this form afterward.


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

@fmorency

This comment has been minimized.

fmorency commented Mar 27, 2015

+1 with @arokem

Would merge the two PR and fine tune on separate PR if needed.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Mar 27, 2015

I don't see where it's broken actually (people just use it in a random fashion), and it might take 2-3 months to be re-merged again if I do another PR. Since it's just a small definition change I'd change it here before forgetting about it.

@arokem

This comment has been minimized.

Member

arokem commented Mar 27, 2015

If it's small, I don't see why it can't be another small PR on top of this.

If you think you might forget to do it, add it to the issues, and assign
yourself to it.

Unless I hear any loud objections, I think that I will merge this tomorrow
or the day after.

On Fri, Mar 27, 2015 at 11:11 AM, Samuel St-Jean notifications@github.com
wrote:

I don't see where it's broken actually (people just use it in a random
fashion), and it might take 2-3 months to be re-merged again if I do
another PR. Since it's just a small definition change I'd change it here
before forgetting about it.


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

@arokem arokem force-pushed the arokem:mdesco-piesno-example branch from bd11490 to 9539072 Mar 29, 2015

arokem added a commit that referenced this pull request Mar 29, 2015

@arokem arokem merged commit ae5a0ff into nipy:master Mar 29, 2015

1 check passed

continuous-integration/travis-ci/pr The Travis CI build passed
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment