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 estimate_sigma bias correction + update example #651

Merged
merged 7 commits into from Jun 5, 2015

Conversation

Projects
None yet
2 participants
@samuelstjean
Contributor

samuelstjean commented May 17, 2015

So I added the correction factor for computing the std of magnitude images. Since they do not include any negative value, the variance is lower and the mean is shifted as compared ot a gaussian distribution. This factor corrects that fact and returns a closer estimate of the variance.
I also updated the example to show this fact and changed it to show a slice of diffusion instea,d since b0 images are always relatively clean.

Since estimate_sigma is the same as computing the variance over the background, the correction factor needs an estimate of the local magnitude SNR, but it is zero in that specific case. Hence, the value can be precomputed without using hyp1f1. Next step is to have it on a voxelwise basis with said hyp1f1 function in need of love over here : #650

samuelstjean added some commits May 17, 2015

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented May 17, 2015

Since I expect the "why did you change the tests" question, an example as to why this way is better than without the correction factor:

% Generate white noise N(25,50**2)
In [34]: a=(np.random.randn(10000)*50)+25

In [35]: np.mean(a)
Out[35]: 24.437325656547163

In [36]: np.std(a)
Out[36]: 50.696824939621159

In [37]: b=(np.random.randn(10000)*50)+25

In [38]: np.std(b)
Out[38]: 49.9460988165076

In [39]: np.mean(b)
Out[39]: 25.756600844271365

% Turn that to Rician noise
In [40]: c=np.sqrt(a**2+b**2)

% Variance is lower, since we lost the negative numbers
In [41]: np.std(c)
Out[41]: 36.449713000167094

% The mean is also shifted, resulting in the signal dependent bias
In [42]: np.mean(c)
Out[42]: 70.68793843681992

%Magic correction factor, which I precomputed with this
In [43]: from scilpy.denoising.stabilizer import _test_xi

% Since we assume zero mean background in this PR, 
% the factor is not signal dependent, but only N coil(s) dependent 
In [44]: _test_xi(25,50,1)
Out[44]: 0.47990987386190764

% Apply the factor, which return a much more closer estimate, 
% (51 instead of 35), while the real value hovers about 50
In [45]: 35.86981271463975/sqrt(0.47990987386190764)
Out[45]: 51.778476330177732
@arokem

This comment has been minimized.

Member

arokem commented May 19, 2015

Cool. I am going to read the Koay paper and give this a more proper review. For now, just a first question: does it make sense to leave the 'old' functionality as an option here?

if N in correction_factor:
factor = correction_factor[N]
else:
factor = correction_factor[1]

This comment has been minimized.

@arokem

arokem May 19, 2015

Member

If the user provides a non-sense value of N, I think that it would make more sense to raise an error, rather than assuming that they meant 1.

This comment has been minimized.

@samuelstjean

samuelstjean May 19, 2015

Contributor

I was thinking of raising a warning, but if you don't force any value all the old code will break in your face.

This comment has been minimized.

@arokem

arokem May 19, 2015

Member

I think that you are forcing a value, by setting the default in the
function signature.

I just meant that if a user called:

estimate_sigma(arr, N=3)

The function should raise an error, not silently behave as though the user
entered N=1. Or am I missing something?

On Tue, May 19, 2015 at 7:57 AM, Samuel St-Jean notifications@github.com
wrote:

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

  • Precomputed factor from Koay 2006, this corrects the bias of magnitude image

  • correction_factor = {1: 0.42920367320510366,
  •                     4: 0.4834941393603609,
    
  •                     6: 0.4891759468548269,
    
  •                     8: 0.49195420135894175,
    
  •                    12: 0.4946862482541263,
    
  •                    16: 0.4960339908122364,
    
  •                    20: 0.4968365823718557,
    
  •                    24: 0.49736907650825657,
    
  •                    32: 0.49803177052530145,
    
  •                    64: 0.49901964176235936}
    
  • if N in correction_factor:
  •    factor = correction_factor[N]
    
  • else:
  •    factor = correction_factor[1]
    

I was thinking of raising a warning, but if you don't force any value all
the old code will break in your face.


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

This comment has been minimized.

@samuelstjean

samuelstjean May 19, 2015

Contributor

No I am indeed forcing a value. I just precomputed a bunch of mostly plausible values. What I'm not sure of though is what to do if the user does not input any value for N. All the old code won't be working anymore, but I could raise an error if the number is not in the list, 1 for no input, and the value as requested if it's in the list.

This comment has been minimized.

@arokem

arokem May 19, 2015

Member

Exactly: raise an error if a number is provided, but it is not a key to the
correction_factor dict. Otherwise, default to 1.

I propose to rewrite this logic as:

if N in correction_factor:
    factor = correction_factor[N]  # Per default, N = 1
else:
   raise ValueError("Correction factor needs to be one of the

following: %s, but you entered %s" %(list(correction_factor.keys()), N))

With the necessary line breaks, etc.

On Tue, May 19, 2015 at 8:37 AM, Samuel St-Jean notifications@github.com
wrote:

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

  • Precomputed factor from Koay 2006, this corrects the bias of magnitude image

  • correction_factor = {1: 0.42920367320510366,
  •                     4: 0.4834941393603609,
    
  •                     6: 0.4891759468548269,
    
  •                     8: 0.49195420135894175,
    
  •                    12: 0.4946862482541263,
    
  •                    16: 0.4960339908122364,
    
  •                    20: 0.4968365823718557,
    
  •                    24: 0.49736907650825657,
    
  •                    32: 0.49803177052530145,
    
  •                    64: 0.49901964176235936}
    
  • if N in correction_factor:
  •    factor = correction_factor[N]
    
  • else:
  •    factor = correction_factor[1]
    

No I am indeed forcing a value. I just precomputed a bunch of mostly
plausible values. What I'm not sure of though is what to do if the user
does not input any value for N. All the old code won't be working anymore,
but I could raise an error if the number is not in the list, 1 for no
input, and the value as requested if it's in the list.


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

This comment has been minimized.

@arokem

arokem May 19, 2015

Member

Just to be clear: I believe that the 'default to 1' part is already taken
care of.

On Tue, May 19, 2015 at 8:44 AM, Ariel Rokem arokem@gmail.com wrote:

Exactly: raise an error if a number is provided, but it is not a key to
the correction_factor dict. Otherwise, default to 1.

I propose to rewrite this logic as:

if N in correction_factor:
    factor = correction_factor[N]  # Per default, N = 1
else:
   raise ValueError("Correction factor needs to be one of the

following: %s, but you entered %s" %(list(correction_factor.keys()), N))

With the necessary line breaks, etc.

On Tue, May 19, 2015 at 8:37 AM, Samuel St-Jean notifications@github.com
wrote:

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

  • Precomputed factor from Koay 2006, this corrects the bias of magnitude image

  • correction_factor = {1: 0.42920367320510366,
  •                     4: 0.4834941393603609,
    
  •                     6: 0.4891759468548269,
    
  •                     8: 0.49195420135894175,
    
  •                    12: 0.4946862482541263,
    
  •                    16: 0.4960339908122364,
    
  •                    20: 0.4968365823718557,
    
  •                    24: 0.49736907650825657,
    
  •                    32: 0.49803177052530145,
    
  •                    64: 0.49901964176235936}
    
  • if N in correction_factor:
  •    factor = correction_factor[N]
    
  • else:
  •    factor = correction_factor[1]
    

No I am indeed forcing a value. I just precomputed a bunch of mostly
plausible values. What I'm not sure of though is what to do if the user
does not input any value for N. All the old code won't be working anymore,
but I could raise an error if the number is not in the list, 1 for no
input, and the value as requested if it's in the list.


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

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented May 19, 2015

Only one would want to compare with an old result, but I take for granted that people don't upgrade libs for fun when doing a huge project. I don't think the old behavior is still needed, since it's only an added division, kind of like when we fixed the nlmeans that was missing a square root.

The real way to do it would be to use the function instead of a precomputed table (that way it will work for all values of N), but that need a working hyp1f1 function, so it's gonna be that for the time being.

@arokem

This comment has been minimized.

Member

arokem commented May 20, 2015

Concerning the main part of this: do I understand correctly that you are simply applying equation 7 of the Koay and Basser paper, for the case where S=0? Might be worth adding to the documentation of the function. Also, how did you compute these actual numbers? The paper only implicitly includes them as the intercepts of the curves in Figure 1.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented May 20, 2015

Yep, that's equation 7, but the factor xi is given by equation 8 for the Rician case and equation 18 for the general case. I indeed used eta=0 (and thus theta=eta/sigma=0, with theta the magnitude SNR, eta the signal value and sigma the noise std). Line 299 depicts this fact with a comment, and the reference is in the function docstring.

-------
This function is the same as manually taking the standard deviation of the
background and gives one value for the whole 3D array. Since it was
developped for T1 imaging. It is expected to perform ok on diffusion MRI data,

This comment has been minimized.

@arokem

arokem May 20, 2015

Member

I think this was meant to be one sentence, and there's a small typo in here ("developped"). It's also a bit vague, so might not help users that much. How about the following:

This function implements equation 7 in [Koay2006]_: it applies a coil-number dependent correction factor to the 
standard deviation of the background and  produces a single value for the input volume, or 4D volume series. Since 
this formulation was developed for T1 imaging, it could work well for diffusion MRI, but take note: when used for 
denoising, this method may oversmooth some regions and leave others un-denoised. This is particularly true for cases 
in which the noise profile is not uniform across the volume. Consider using :func:`piesno` if visual examination  reveals 
non-uniform denoising performance.

As an aside, the Notes section should precede the "References" section.

@@ -292,6 +331,6 @@ def estimate_sigma(arr, disable_background_masking=False):
for i in range(sigma.size):
convolve(arr[..., i], k, output=conv_out)
mean_block = np.sqrt(6/7) * (arr[..., i] - 1/6 * conv_out)
sigma[i] = np.sqrt(np.mean(mean_block[mask]**2))
sigma[i] = np.sqrt(np.mean(mean_block[mask]**2) / factor)

This comment has been minimized.

@arokem

arokem May 20, 2015

Member

PEP8: white-space around operators (I know - it was there before, but let's please clean it up, while we're at it).

This comment has been minimized.

@samuelstjean

samuelstjean May 20, 2015

Contributor

Wasn't there a note saying you can keep the priority of operation grouped together like in

a**2 + b**2

Plus I really find it counterintuitive to have the exponent not close to the actual number, like in

a ** 2 / b ** 2

Seems messy :/

This comment has been minimized.

@arokem

arokem May 20, 2015

Member

Yes: operators with higher precedence should be without white-space. But I
don't think this applies in cases where they are parenthesis-delimited,
such as here.

All the BDFL gave us mere mortals is the following guidance:

Yes:

i = i + 1
submitted += 1
x = x_2 - 1
hypot2 = x_x + y*y
c = (a+b) * (a-b)

No:

i=i+1
submitted +=1
x = x * 2 - 1
hypot2 = x * x + y * y

c = (a + b) * (a - b)

I don't know what to say, except that an automatic tool (like pep8) would
definitely flag this one. I agree with you that it doesn't read as nicely,
but that's the good thing about using standards: no one's individual taste
matters.

On Tue, May 19, 2015 at 6:03 PM, Samuel St-Jean notifications@github.com
wrote:

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

@@ -292,6 +331,6 @@ def estimate_sigma(arr, disable_background_masking=False):
for i in range(sigma.size):
convolve(arr[..., i], k, output=conv_out)
mean_block = np.sqrt(6/7) * (arr[..., i] - 1/6 * conv_out)

  •    sigma[i] = np.sqrt(np.mean(mean_block[mask]**2))
    
  •    sigma[i] = np.sqrt(np.mean(mean_block[mask]**2) / factor)
    

Wasn't there a note saying you can keep the priority of operation grouped
together like in

a_2 + b_2

Plus I really find it counterintuitive to have the exponent not close to
the actual number, like in

a ** 2 / b ** 2
Seems messy :/


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

This comment has been minimized.

@samuelstjean

samuelstjean May 20, 2015

Contributor

PEP 8:

But most importantly: know when to be inconsistent -- sometimes the style guide just doesn't apply. When in doubt, use your best judgment.

This comment has been minimized.

@arokem

arokem May 20, 2015

Member

And if it wasn't bad enough, he threw that one in for good measure. I
wonder how many programmer-hours have gone to debating that one line over
the years
On May 19, 2015 6:24 PM, "Samuel St-Jean" notifications@github.com wrote:

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

@@ -292,6 +331,6 @@ def estimate_sigma(arr, disable_background_masking=False):
for i in range(sigma.size):
convolve(arr[..., i], k, output=conv_out)
mean_block = np.sqrt(6/7) * (arr[..., i] - 1/6 * conv_out)

  •    sigma[i] = np.sqrt(np.mean(mean_block[mask]**2))
    
  •    sigma[i] = np.sqrt(np.mean(mean_block[mask]**2) / factor)
    

PEP 8:

But most importantly: know when to be inconsistent -- sometimes the style
guide just doesn't apply. When in doubt, use your best judgment.


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

This comment has been minimized.

@arokem

arokem May 20, 2015

Member

But enough fun and games. Let me know when you've addressed the comments
and I'll take another look
On May 19, 2015 6:28 PM, "Ariel Rokem" arokem@gmail.com wrote:

And if it wasn't bad enough, he threw that one in for good measure. I
wonder how many programmer-hours have gone to debating that one line over
the years
On May 19, 2015 6:24 PM, "Samuel St-Jean" notifications@github.com
wrote:

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

@@ -292,6 +331,6 @@ def estimate_sigma(arr, disable_background_masking=False):
for i in range(sigma.size):
convolve(arr[..., i], k, output=conv_out)
mean_block = np.sqrt(6/7) * (arr[..., i] - 1/6 * conv_out)

  •    sigma[i] = np.sqrt(np.mean(mean_block[mask]**2))
    
  •    sigma[i] = np.sqrt(np.mean(mean_block[mask]**2) / factor)
    

PEP 8:

But most importantly: know when to be inconsistent -- sometimes the style
guide just doesn't apply. When in doubt, use your best judgment.


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

arr = np.zeros((3, 3, 3))
arr[0, 0, 0] = 1
sigma = estimate_sigma(arr, disable_background_masking=True)
assert_array_almost_equal(sigma, 0.46291005)
sigma = estimate_sigma(arr, disable_background_masking=True, N=4)

This comment has been minimized.

@arokem

arokem May 20, 2015

Member

Nice. Good to test for numbers other than 1...

@arokem

This comment has been minimized.

Member

arokem commented May 20, 2015

Any ideas how this would affect the use of estimate_sigma in the RESTORE example?

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented May 20, 2015

It's almost the same as using the multiplication by sqrt(pi/2) I'd say, but this one is more general since it takes into account the increased noise floor with higher N. Good thing you brought that up, I think it's better to untangle the noise estimation method from the models, else two correction factors would get applied.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented May 20, 2015

Fixed in majority

@arokem

This comment has been minimized.

arokem commented on dipy/denoise/noise_estimate.py in 3b5f9ae May 20, 2015

[1] => [1]_

@arokem

This comment has been minimized.

arokem commented on dipy/denoise/noise_estimate.py in 3b5f9ae May 20, 2015

[2] => [2]_

@arokem

This comment has been minimized.

arokem commented on dipy/denoise/noise_estimate.py in 3b5f9ae May 20, 2015

Please format the references so that they get linked to the rest above. See for example here the top-level docstring here: https://github.com/nipy/dipy/blob/master/dipy/reconst/sfm.py#L2 for the syntax.

@arokem

This comment has been minimized.

Member

arokem commented on 9d40363 May 20, 2015

🆒

@arokem

This comment has been minimized.

Member

arokem commented May 20, 2015

OK - I think that we can merge this. Let's plan on some time next week, so that others can take a look as well. Anyone else have any comments on this?

@arokem

This comment has been minimized.

Member

arokem commented May 20, 2015

Sorry - I just remembered one more question that I had: what effect does this have on the RESTORE example? Did you run that one?

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented May 20, 2015

nope, but I can absolutely garantee that we should remove the sqrt(pi/2) fix, since that overlaps with this correction factor. Hence the now deeply buried issue I posted with my suggested fix (this is number 1 of them if I recall correctly).

@arokem

This comment has been minimized.

Member

arokem commented May 20, 2015

OK - since this will have an effect there, please do run it and make sure that it runs, make any corrections to that example that you think are necessary, as part of this PR. Which one is the 'deeply buried' issue you refer to (link)?

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented May 20, 2015

@arokem

This comment has been minimized.

Member

arokem commented May 20, 2015

Oh - that issue - all good suggestions, as discussed. So you think the
first on that list should be a part of this PR? Seems like a good idea to
me. Either way, we can't merge this if it breaks the RESTORE example.

On Wed, May 20, 2015 at 10:14 AM, Samuel St-Jean notifications@github.com
wrote:

#636 #636


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

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented May 20, 2015

So, do you know the specs of the scanner used for the standford dataset?

@arokem

This comment has been minimized.

Member

arokem commented May 20, 2015

Yes:

http://cni.stanford.edu/wiki/MR_Scanner

We used the 32 channel coil for these measurements.

On Wed, May 20, 2015 at 10:57 AM, Samuel St-Jean notifications@github.com
wrote:

So, do you know the specs of the scanner used for the standford dataset?


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

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented May 20, 2015

Well, apparently GE can use both ARC (GRAPPA) or SENSE from a quick internet search. Any idea which one is used by any chance?

@arokem

This comment has been minimized.

Member

arokem commented May 20, 2015

99% sure it was SENSE

On Wed, May 20, 2015 at 11:13 AM, Samuel St-Jean notifications@github.com
wrote:

Well, apparently GE can use both ARC (GRAPPA) or SENSE from a quick
internet search. Any idea which one is used by any chance?


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

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented May 20, 2015

I took a look at the code, and apparently everything is fine, I thought the noise estimation was embedded in the method, but you use another estimator, which is fine and completely independent from this implementation. I orignally though two different correction factor would be applied, but that's not the case.

@arokem

This comment has been minimized.

Member

arokem commented May 20, 2015

Thanks for all the cleanup in dti.py!

This function is being used in the RESTORE example:

https://github.com/nipy/dipy/blob/master/doc/examples/restore_dti.py#L162

Could you please check what happens when you run this example?

@arokem

This comment has been minimized.

Member

arokem commented May 22, 2015

Just checking back in here: did you have a chance to run the RESTORE example on your branch?

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented May 22, 2015

Not yet, I just learned I had to defend my master thesis by next week (which I finished wednesday), yesterday. So that will probably take my time until then.

@arokem

This comment has been minimized.

Member

arokem commented May 22, 2015

Oh wow! This can definitely wait... Hopefully more people will take a look
in the meanwhile.

Good luck!

On Fri, May 22, 2015 at 10:39 AM, Samuel St-Jean notifications@github.com
wrote:

Not yet, I just learned I had to defend my master thesis by next week
(which I finished wednesday), yesterday. So that will probably take my time
until then.


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

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented May 27, 2015

Just ran it, doesn't seem to change anything visually (as expected). Fuunnily enough, the fmri guys on our side hardcoded a bunch of values since they were unhappy with this function result, and with the proper correction it gives the same value that they hardcoded, for whatever that is worth.
Well, maybe anyone else could comment on the usefulness of changing this as default, but as far as I am concerned this is the right thing to do. And it's also the first open source implementation that I am aware of ;)

dti_fa_distributions
tensor_ellipsoids_restore_noisy
tensor_ellipsoids_wls
tensor_ellipsoids_wls_noisy

@arokem

This comment has been minimized.

Member

arokem commented May 28, 2015

Well - if this will make your fMRI people happy, then I am happy too! :-)

Thanks for checking that. This is complete from my point-of-view, but I
would like to give others a chance to take a look. If no-one pipes up, I
will merge this next week.

On Wed, May 27, 2015 at 4:04 PM, Samuel St-Jean notifications@github.com
wrote:

Just ran it, doesn't seem to change anything visually (as expected).
Fuunnily enough, the fmri guys on our side hardcoded a bunch of values
since they were unhappy with this function result, and with the proper
correction it gives the same value that they hardcoded, for whatever that
is worth.
Well, maybe anyone else could comment on the usefulness of changing this
as default, but as far as I am concerned this is the right thing to do. And
it's also the first open source implementation that I am aware of ;)

[image: dti_fa_distributions]
https://cloud.githubusercontent.com/assets/3030760/7849388/dbe4d568-04a2-11e5-9bd3-4684beab4f03.png
[image: tensor_ellipsoids_restore_noisy]
https://cloud.githubusercontent.com/assets/3030760/7849391/dbeab168-04a2-11e5-8b27-8475aebad7d7.png
[image: tensor_ellipsoids_wls]
https://cloud.githubusercontent.com/assets/3030760/7849390/dbe9fae8-04a2-11e5-88c6-ab3d2a17fe9c.png
[image: tensor_ellipsoids_wls_noisy]
https://cloud.githubusercontent.com/assets/3030760/7849389/dbe9938c-04a2-11e5-86f2-28c6808dcda4.png


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

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented May 28, 2015

This broke for an out of space error on travis, no need to retrigger it.

s1
s2

@arokem

This comment has been minimized.

Member

arokem commented May 28, 2015

Interesting. Haven't seen that one before...

On Wed, May 27, 2015 at 8:22 PM, Samuel St-Jean notifications@github.com
wrote:

This broke for an out of space error on travis, no need to retrigger it.

[image: s1]
https://cloud.githubusercontent.com/assets/3030760/7852163/30a86104-04c7-11e5-9a82-cd7a1312f53b.png
[image: s2]
https://cloud.githubusercontent.com/assets/3030760/7852164/30b0c272-04c7-11e5-9ca3-c779e7b583e0.png


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

arokem added a commit that referenced this pull request Jun 5, 2015

Merge pull request #651 from samuelstjean/nlmeans_bias_correct
Added estimate_sigma bias correction + update example

@arokem arokem merged commit 2aed9ea into nipy:master Jun 5, 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