[MRG+1] more closesly match the BayesShrink paper in _wavelet_threshold

[grlee77]

Description

pinging @stsievert, @JDWarner, @sciunto, @soupault from the discussion in #2190 for feedback on this one.

This small PR modifies the recently merged #2190 to more closely match the implementation of the Chang et. al. IEEE paper mentioned in the docstring. I apologize for not being more active in reviewing that PR before it was merged. I think either we need to make these changes to match the reference, or we need to change the docstring to match what is currently implemented.

I realize there is also ongoing work in #2240 related to color denoising.

There are no API changes proposed, but the resulting denoised images look much different than before.

The primary changes made to match the referenced article are:

• approximation coefficients are not soft thresholded
• BayesShrink thresholds are computed individually for each individual detail coefficient sub-band
• the wavelet decomposition is not performed to the full possible number of levels (although I did not fix at 4 levels as in the reference as we should support any size of input images)

I have not tested extensively whether using separate thresholds per sub-band is truly better, but if we cite the BayesShrink paper in the References we should more closely match the implementation in the citation.

Demonstration

Here is a quick demo for which I have pasted the image for current master as opposed to this PR.

import numpy as np
from matplotlib import pyplot as plt

from skimage import img_as_float
import skimage.data
from skimage.restoration._denoise import _wavelet_threshold

x = img_as_float(skimage.data.camera())
sigma = 0.1
y = x + sigma*np.random.randn(*x.shape)
ydn = _wavelet_threshold(y, 'db1')

fig, axes = plt.subplots(1, 2)
axes[0].imshow(y, interpolation='nearest')
axes[0].set_title('noisy')
axes[0].set_xticks([])
axes[0].set_yticks([])
axes[1].imshow(ydn, interpolation='nearest')
axes[1].set_title('with proposed changes')
axes[1].set_xticks([])
axes[1].set_yticks([])
plt.tight_layout()

Result with code in this PR:

or when providing sigma=0.9*sigma as anargument to `_wavelet_threshold:

ydn_fixed_thresh = _wavelet_threshold(y, 'db1', sigma=0.9*sigma)

It is still possible to supply threshold to get a single threshold used at all levels:

ydn_fixed_thresh = _wavelet_threshold(y, 'db1', threshold=sigma)

The default level of denoising produced by the master branch is much lower:

[codecov-io]

Current coverage is 90.63% (diff: 96.92%)

Merging #2241 into master will increase coverage by 0.02%

@@             master      #2241   diff @@
==========================================
  Files           304        304          
  Lines         21475      21521    +46   
  Methods           0          0          
  Messages          0          0          
  Branches       1846       1850     +4   
==========================================
+ Hits          19458      19505    +47   
  Misses         1663       1663          
+ Partials        354        353     -1   

Powered by Codecov. Last update 92a3851...c4d8058

[soupault]

Please, move this to the docstring Notes section.

[soupault]

@grlee77 thanks for the work done! I have just a couple of minor stylistic comments.
I personally prefer to merge this after #2242, and it might require some extra work.

[grlee77]

@soupault. I agree with the suggestions. waiting for #2242 sounds reasonable.

[grlee77]

travis failure was just a time-out on OS X

[sciunto]

@grlee77 Could you rebase please?

[soupault]

Shouldn't we expose this parameter in the outmost function?

[soupault]

I've made some remarks, but this PR looks good to me in general.

[sciunto]

@grlee77 kind ping. Do you have a chance to look at @soupault comments? I think it's closed to be merged. ;)

[grlee77]

@sciunto, @soupault
I had to rebase again, but the last 3 commits are the new ones. I think they address the remaining comments.

[soupault]

@grlee77 Please, set multichannel to be the last argument.

[soupault]

One more comment, then 👍.
Thanks @grlee77!

[grlee77]

done. I also avoided duplicate code for the mean-absolute-deviation estimate of sigma by reusing _sigma_est_dwt (I moved the _sigma_est_dwt function's location above _wavelet_threshold so that it is defined prior to its first use).

Is there a policy/convention that multichannel should always be the last argument? If so, estimate_sigma and denoise_bilateral are not currently following that convention, but I didn't change them here.

Also, on another tangent, I noticed that all of the denoising tests use assert rather than the assert_ function from numpy.testing. IIRC assert calls can be skipped if python is called with -O so I thought assert_ was generally preferred. I can make a small separate PR for that if it is of interest.

[soupault]

@grlee77 Perfect, thanks!
No, I don't think there is an actual convention on that matter. Just from my perspective, there are 3 general keywords that we may finish adding to the majority of public functions: multichannel, preserve_range and out. It would be better to have them as the last parameters from the start. As dwt is a new feature and hasn't yet been released (so the API isn't stabilized yet), I think, you could apply this rule to other functions as well.
For denoise_bilateral, I suppose, this is not so pleasant due to the deprecation stuff, so let's just leave it there as is.
P.S. Don't forget to change the order in the docstring as well. 😉.

Regarding asserts: yes from me to change.

[jni]

@grlee77 thank you!