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 gaussian noise option to estimate_sigma #699

Merged
merged 4 commits into from Oct 7, 2015

Conversation

Projects
None yet
5 participants
@samuelstjean
Contributor

samuelstjean commented Aug 11, 2015

Since this is a case that can happen, and was the old behavior, I added an option to disable the correction factor, which is the default as it was before.

@arokem

This comment has been minimized.

arokem commented on dipy/denoise/noise_estimate.py in 5a9d314 Aug 11, 2015

Invalid syntax. You need a comma after 1

@arokem

This comment has been minimized.

Member

arokem commented Aug 11, 2015

Would be good to have a test that checks this mode.

samuelstjean added some commits Aug 11, 2015

@arokem

This comment has been minimized.

Member

arokem commented Aug 17, 2015

LGTM. Anyone else have any comments here? If not, I will merge this in a few days.

@arokem

This comment has been minimized.

Member

arokem commented Aug 25, 2015

Hey @mdesco - any opinions here? I don't know this stuff well enough, and could really use a second opinion.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Aug 27, 2015

For future reference, this merely restores the default behavior, but also keeps the suggested enhancements for the fmri people.

@mdesco

This comment has been minimized.

Contributor

mdesco commented Aug 27, 2015

Not much to say. I asked @samuelstjean to put the N=0 option for the Gaussian noise case. The MGH group now has a paper to correct the phase and produce DWI with Gaussian-like profiles. Hence, the N=0 case can become important to have.

http://www.ncbi.nlm.nih.gov/pubmed/26241680

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Aug 27, 2015

@mdesco should that be the default one?

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Aug 27, 2015

Or N=1 is more likely?

@mdesco

This comment has been minimized.

Contributor

mdesco commented Aug 27, 2015

To me, N=1 should be the default one, the rician case.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Aug 27, 2015

Depends, if you have fancy data, N=1 is probably advisable, N=0 is the previous default, and the one used in the paper. So to keep the original behavior, that is also the default. Since this is a data driven estimation, this keeps the spirit of the original version (well, actually it should do exactly as the original provided we did not introduce a bug or anything in the rewrite).

Remember this is for the estimation of the standard deviation (with a n additional correction factor as needed), while the rician correction in nlmeans is set to true by default)

Number of coils of the receiver array. Use N = 1 in case of a SENSE
reconstruction (Philips scanners) or the number of coils for a GRAPPA
reconstruction (Siemens and GE). See [1] for more information.
reconstruction (Siemens and GE). Use 0 to disable the correction factor,
as for example is the noise is Gaussian distributed. See [1] for more information.

This comment has been minimized.

@arokem

arokem Aug 28, 2015

Member

=> "Use 0 to disable the correction factor, for example if the noise is Gaussian distributed"

@arokem

This comment has been minimized.

Member

arokem commented Aug 28, 2015

Hi @samuelstjean - what do you mean by "fancy"? 2006's "fancy" is today's extreme low end.

If by "fancy" you mean a multi-coil or accelerated acquisition than I am with @mdesco on this: 1 should be the default.

@rfdougherty

This comment has been minimized.

rfdougherty commented Aug 28, 2015

I agree with default=1. I don't think any standard vendor sequence produces real data with Gaussian noise.

Also, note that the current EPI diffusion sequences on GE use a SENSE-like recon (they call it ASSET) for inplane acceleration. (Siemens is currently the only one of the three that uses GRAPPA by default for inplane acceleration.) On GE, the appropriate N is 1 if ASSET (i.e. SENSE) is used, or num coils if ASSET is disabled (with ASSET disabled, coils are combined with a simple sum-of-squares algorithm). I suspect Philips would work the same (but can't say for sure).

So the recommended settings need to be more nuanced: "For stock vendor sequences, you probably want to use N=num coils on Siemens. On GE and Philips, you probably want N=1 if inplane acceleration is enabled, otherwise use N=num coils. But it's best to confirm the appropriate N to use by inspecting your actual noise distribution." (That disclaimer is essential, since any of this can change. E.g., GE is saying they'll release a GRAPPA-like option for dw-epi soon, Siemens may add SENSE coil combination, etc.)

And maybe even add "If you use non-stock sequences, you should talk to the pulse sequence developers for advice." E.g., at our center, we do GRAPPA followed by SENSE coil combination (a.k.a. "SENSE1"), and phase correction to produce real data with zero-mean gaussian noise. So for us, N=0 would be appropriate, even though we use a GE scanner. (And as @mdesco notes, this will be more common given the recent paper out of MGH.)

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Aug 28, 2015

That seems much more confusing than needed for anyone having no idea of what they are supposed to use. Plus, even if you ask for the appropriate answer, the number people who can help you in any given institute which might be able to help you can be counted on one hand (if they are cooperative and do exist in the first place).

I still vote for restoring the default option as it was before (N=0), and people needing it can still change the appropriate value from N=0 to whatever. Remember I'm only trying to shove something not belonging there in the original design for more accurate estimation, and people complained from that, so hence the "let's return to the default value".

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Aug 28, 2015

By fancy I meant anything that deviates from the standard, stationary Rician/non central chi noise model, so your everyday acquisition would be a not fancy one while the HCP/multiband/etc. and the just released paper on phase smoothing for gaussian distributed noise is considered fancy in my book.

I think that people are confusing two un-related but still close issues here : the rician intensity correction bias from nlmeans and the second order moment correction factor (or variance correction for non gaussian data) from estimate_sigma.

This PR is about the second case, that is the old estimate_sigma estimator is a statistical estimator designed for gaussian noise and I tacked together a correction factor to account for the lower variance arising from not having a gaussian distribution. Of course this is not what the original version, so I just put back the original estimation instead of using that correction factor by default as orginally intended. This is also what the current release does.

Where I think the mix-up arises from, and that is totally untouched and still intact, is that nlmeans use the rician bias correction factor by default by substracting -2*variance to each voxel at the far end of the denoising process to correct the bias intensity. This is still the case and this is also what I think people wrongly refer to as N=1 should be default. You can look at the ~7 lines I changed to convince yourself of the good natured idea of the proposed change of course.

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Sep 14, 2015

Okay @rfdougherty and @mdesco and all. I think Sam gave a complete explanation about this. If you do not disagree I would like to go ahead and merge this PR. Please let me know otherwise.

As a future project it would be nice to be able to detect these distributions automatically. But this is good for now.

You have 48 hours for asking for extra time to look into this PR. If nobody responds I will go ahead and merge.

@rfdougherty

This comment has been minimized.

rfdougherty commented Sep 14, 2015

OK with me.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Sep 21, 2015

As for detecting the noise distribution, I have an idea on how to do that, maybe, I will need to play a bit with that beforehand, but it should be doable with some maths.

@arokem

This comment has been minimized.

Member

arokem commented Oct 2, 2015

@Garyfallidis - did you mean to merge this one?

Garyfallidis added a commit that referenced this pull request Oct 7, 2015

Merge pull request #699 from samuelstjean/patch-2
Added gaussian noise option to estimate_sigma

@Garyfallidis Garyfallidis merged commit 1cabdb4 into nipy:master Oct 7, 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