From d41f465ce5ba9d7eb22122ea3125697fa514465f Mon Sep 17 00:00:00 2001
From: "Gregory R. Lee" <gregory.lee@cchmc.org>
Date: Wed, 10 Aug 2016 15:57:29 -0400
Subject: [PATCH 01/16] DOC: correct SUREShrink->BayesShrink in the
 _wavelet_threshold docstring

---
 skimage/restoration/_denoise.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index 0e9f6df8bcd..15f51df17f8 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -353,8 +353,8 @@ def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft'):
         is None (the default) by the method in [2]_.
     threshold : float, optional
         The thresholding value. All wavelet coefficients less than this value
-        are set to 0. The default value (None) uses the SureShrink method found
-        in [1]_ to remove noise.
+        are set to 0. The default value (None) uses the BayesShrink method
+        found in [1]_ to remove noise.
     mode : {'soft', 'hard'}, optional
         An optional argument to choose the type of denoising performed. It
         noted that choosing soft thresholding given additive noise finds the

From 5d9dba3e1cd4c6b44a89028234d7377dd028de47 Mon Sep 17 00:00:00 2001
From: "Gregory R. Lee" <gregory.lee@cchmc.org>
Date: Wed, 10 Aug 2016 15:58:23 -0400
Subject: [PATCH 02/16] FIX: _wavelet_threshold: do not denoise the
 approximation coefficients

---
 skimage/restoration/_denoise.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index 15f51df17f8..099ce4a7d0f 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -389,8 +389,7 @@ def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft'):
 
     denoised_detail = [{key: pywt.threshold(level[key], value=threshold,
                        mode=mode) for key in level} for level in coeffs[1:]]
-    denoised_root = pywt.threshold(coeffs[0], value=threshold, mode=mode)
-    denoised_coeffs = [denoised_root] + [d for d in denoised_detail]
+    denoised_coeffs = [coeffs[0]] + [d for d in denoised_detail]
     return pywt.waverecn(denoised_coeffs, wavelet)
 
 

From 799a21efa211907c19d3afdd611ac7d70774f462 Mon Sep 17 00:00:00 2001
From: "Gregory R. Lee" <gregory.lee@cchmc.org>
Date: Wed, 10 Aug 2016 16:11:13 -0400
Subject: [PATCH 03/16] ENH: BayesShrink thresholds are computed for each
 coefficient array

This matches what is done in the Chang et. al. reference from the DocString.
---
 skimage/restoration/_denoise.py | 33 +++++++++++++++++++++++++++------
 1 file changed, 27 insertions(+), 6 deletions(-)

diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index 099ce4a7d0f..5936f8007e1 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -335,6 +335,15 @@ def denoise_tv_chambolle(im, weight=0.1, eps=2.e-4, n_iter_max=200,
     return out
 
 
+def _bayes_thresh(details, var):
+    """BayesShrink threshold for a zero-mean details coeff array."""
+    # equivalent to:  dvar = np.var(details) for 0-mean details array
+    dvar = np.mean(details*details)
+    eps = np.finfo(details.dtype).eps
+    thresh = var / np.sqrt(max(dvar - var, eps))
+    return thresh
+
+
 def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft'):
     """Performs wavelet denoising.
 
@@ -384,12 +393,24 @@ def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft'):
         sigma = np.median(np.abs(detail_coeffs)) / 0.67448975019608171
 
     if threshold is None:
-        # The BayesShrink threshold from [1]_ in docstring
-        threshold = sigma**2 / np.sqrt(max(img.var() - sigma**2, 0))
-
-    denoised_detail = [{key: pywt.threshold(level[key], value=threshold,
-                       mode=mode) for key in level} for level in coeffs[1:]]
-    denoised_coeffs = [coeffs[0]] + [d for d in denoised_detail]
+        # The BayesShrink thresholds from [1]_ in docstring
+        var = sigma**2
+        threshold = [{key: _bayes_thresh(level[key], var) for key in level}
+                     for level in coeffs[1:]]
+
+    if np.isscalar(threshold):
+        # a single threshold for all coefficient arrays
+        denoised_detail = [{key: pywt.threshold(level[key],
+                                                value=threshold,
+                                                mode=mode) for key in level}
+                           for level in coeffs[-1]]
+    else:
+        # dict of unique threshold coefficients for each detail coeff. array
+        denoised_detail = [{key: pywt.threshold(level[key],
+                                                value=thresh[key],
+                                                mode=mode) for key in level}
+                           for thresh, level in zip(threshold, coeffs[-1])]
+    denoised_coeffs = [coeffs[0]] + denoised_detail
     return pywt.waverecn(denoised_coeffs, wavelet)
 
 

From 99a1ce8750713dc961c1999afc66738316d7277c Mon Sep 17 00:00:00 2001
From: "Gregory R. Lee" <gregory.lee@cchmc.org>
Date: Wed, 10 Aug 2016 16:17:02 -0400
Subject: [PATCH 04/16] ENH: avoid BayesShrink thresholding on extremely coarse
 resolution coefficient arrays.

This is necessary because these arrays may be extremely small leading to unreliable variance estimation and
noticeable blocklike artifacts in the denoised image.
---
 skimage/restoration/_denoise.py | 23 ++++++++++++++++++-----
 1 file changed, 18 insertions(+), 5 deletions(-)

diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index 5936f8007e1..b9e96f8c685 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -385,31 +385,44 @@ def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft'):
            DOI: 10.1093/biomet/81.3.425
 
     """
-    coeffs = pywt.wavedecn(img, wavelet=wavelet)
-    detail_coeffs = coeffs[-1]['d' * img.ndim]
+    if not isinstance(wavelet, pywt.Wavelet):
+        wavelet = pywt.Wavelet(wavelet)
+
+    # determine the number of wavelet decomposition levels
+    dlen = wavelet.dec_len
+    level = np.min([pywt.dwt_max_level(s, dlen) for s in img.shape])
+    if threshold is None:
+        # Bayes shrink doesn't work well on levels with extremely small
+        # coefficient arrays, so skip a couple of the coarsest levels
+        level = max(level - 2, 1)
+
+    coeffs = pywt.wavedecn(img, wavelet=wavelet, level=level)
+    # detail coefficients at each decomposition level
+    dcoeffs = coeffs[1:]
 
     if sigma is None:
         # Estimates via the noise via method in [2]
+        detail_coeffs = dcoeffs[-1]['d' * img.ndim]
         sigma = np.median(np.abs(detail_coeffs)) / 0.67448975019608171
 
     if threshold is None:
         # The BayesShrink thresholds from [1]_ in docstring
         var = sigma**2
         threshold = [{key: _bayes_thresh(level[key], var) for key in level}
-                     for level in coeffs[1:]]
+                     for level in dcoeffs]
 
     if np.isscalar(threshold):
         # a single threshold for all coefficient arrays
         denoised_detail = [{key: pywt.threshold(level[key],
                                                 value=threshold,
                                                 mode=mode) for key in level}
-                           for level in coeffs[-1]]
+                           for level in dcoeffs]
     else:
         # dict of unique threshold coefficients for each detail coeff. array
         denoised_detail = [{key: pywt.threshold(level[key],
                                                 value=thresh[key],
                                                 mode=mode) for key in level}
-                           for thresh, level in zip(threshold, coeffs[-1])]
+                           for thresh, level in zip(threshold, dcoeffs)]
     denoised_coeffs = [coeffs[0]] + denoised_detail
     return pywt.waverecn(denoised_coeffs, wavelet)
 

From 26067273ef114a1ce067afbac3f82a176e3742f2 Mon Sep 17 00:00:00 2001
From: "Gregory R. Lee" <gregory.lee@cchmc.org>
Date: Wed, 10 Aug 2016 17:38:11 -0400
Subject: [PATCH 05/16] ENH: make wavelet_levels a parameter to
 _wavelet_threshold

---
 skimage/restoration/_denoise.py | 28 +++++++++++++++++++---------
 1 file changed, 19 insertions(+), 9 deletions(-)

diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index b9e96f8c685..25544359834 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -344,7 +344,8 @@ def _bayes_thresh(details, var):
     return thresh
 
 
-def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft'):
+def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft',
+                       wavelet_levels=None):
     """Performs wavelet denoising.
 
     Parameters
@@ -368,6 +369,10 @@ def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft'):
         An optional argument to choose the type of denoising performed. It
         noted that choosing soft thresholding given additive noise finds the
         best approximation of the original image.
+    wavelet_levels : int or None, optional
+        The number of wavelet decomposition levels to use.  The default is
+        three less than the maximum number of possible decomposition levels
+        (where the max depends on the size of `img` and the wavelet to use).
 
     Returns
     -------
@@ -389,14 +394,19 @@ def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft'):
         wavelet = pywt.Wavelet(wavelet)
 
     # determine the number of wavelet decomposition levels
-    dlen = wavelet.dec_len
-    level = np.min([pywt.dwt_max_level(s, dlen) for s in img.shape])
-    if threshold is None:
-        # Bayes shrink doesn't work well on levels with extremely small
-        # coefficient arrays, so skip a couple of the coarsest levels
-        level = max(level - 2, 1)
-
-    coeffs = pywt.wavedecn(img, wavelet=wavelet, level=level)
+    if wavelet_levels is None:
+        # determine the maximum number of possible levels for img
+        dlen = wavelet.dec_len
+        wavelet_levels = np.min(
+            [pywt.dwt_max_level(s, dlen) for s in img.shape])
+
+        # BayesShrink variance estimation doesn't work well on levels with
+        # extremely small coefficient arrays, so skip a few of the coarsest
+        # levels.
+        # Note: ref [1] used a fixed wavelet_levels = 4
+        wavelet_levels = max(wavelet_levels - 3, 1)
+
+    coeffs = pywt.wavedecn(img, wavelet=wavelet, level=wavelet_levels)
     # detail coefficients at each decomposition level
     dcoeffs = coeffs[1:]
 

From d866c999997fd05f66acd69811ab3b8d11623ac9 Mon Sep 17 00:00:00 2001
From: "Gregory R. Lee" <gregory.lee@cchmc.org>
Date: Wed, 10 Aug 2016 18:59:54 -0400
Subject: [PATCH 06/16] TST: fix test:  energy of x is np.sum(x**2) not
 x.sum()**2

---
 skimage/restoration/tests/test_denoise.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py
index 2f193494d82..6135e6828a9 100644
--- a/skimage/restoration/tests/test_denoise.py
+++ b/skimage/restoration/tests/test_denoise.py
@@ -311,9 +311,10 @@ def test_no_denoising_for_small_h():
 
 
 def test_wavelet_denoising():
+    rstate = np.random.RandomState(1234)
     for img, multichannel in [(astro_gray, False), (astro, True)]:
         sigma = 0.1
-        noisy = img + sigma * np.random.randn(*(img.shape))
+        noisy = img + sigma * rstate.randn(*(img.shape))
         noisy = np.clip(noisy, 0, 1)
 
         # Verify that SNR is improved when true sigma is used
@@ -335,17 +336,18 @@ def test_wavelet_denoising():
                                            multichannel=multichannel)
         res2 = restoration.denoise_wavelet(noisy, sigma=sigma,
                                            multichannel=multichannel)
-        assert (res1.sum()**2 <= res2.sum()**2)
+        assert np.sum(res1**2) <= np.sum(res2**2)
 
 
 def test_wavelet_denoising_nd():
+    rstate = np.random.RandomState(1234)
     for ndim in range(1, 5):
         # Generate a very simple test image
         img = 0.2*np.ones((16, )*ndim)
         img[[slice(5, 13), ] * ndim] = 0.8
 
         sigma = 0.1
-        noisy = img + sigma * np.random.randn(*(img.shape))
+        noisy = img + sigma * rstate.randn(*(img.shape))
         noisy = np.clip(noisy, 0, 1)
 
         # Verify that SNR is improved with internally estimated sigma

From 6c512a091ef1e51c86520f60ee12b81fa4fcb9b0 Mon Sep 17 00:00:00 2001
From: "Gregory R. Lee" <gregory.lee@cchmc.org>
Date: Thu, 11 Aug 2016 10:41:20 -0400
Subject: [PATCH 07/16] STY: update style to common skimage conventions

---
 skimage/restoration/_denoise.py | 33 ++++++++++++++++++++-------------
 1 file changed, 20 insertions(+), 13 deletions(-)

diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index 25544359834..2f9f9fbcc08 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -337,7 +337,7 @@ def denoise_tv_chambolle(im, weight=0.1, eps=2.e-4, n_iter_max=200,
 
 def _bayes_thresh(details, var):
     """BayesShrink threshold for a zero-mean details coeff array."""
-    # equivalent to:  dvar = np.var(details) for 0-mean details array
+    # Equivalent to:  dvar = np.var(details) for 0-mean details array
     dvar = np.mean(details*details)
     eps = np.finfo(details.dtype).eps
     thresh = var / np.sqrt(max(dvar - var, eps))
@@ -372,13 +372,24 @@ def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft',
     wavelet_levels : int or None, optional
         The number of wavelet decomposition levels to use.  The default is
         three less than the maximum number of possible decomposition levels
-        (where the max depends on the size of `img` and the wavelet to use).
+        (see Notes below).
 
     Returns
     -------
     out : ndarray
         Denoised image.
 
+    Notes
+    -----
+    Reference [1]_ used four levels of wavelet decomposition.  To be more
+    flexible for a range of input sizes, the implementation here stops 3 levels
+    prior to the maximum level of decomposition for `img` (the exact # of
+    levels thus depends on `img.shape` and the chosen wavelet). BayesShrink
+    variance estimation doesn't work well on levels with extremely small
+    coefficient arrays.  This is the rationale for skipping a few of the
+    coarsest levels.  The user can override the automated setting by explicitly
+    specifying `wavelet_levels`.
+
     References
     ----------
     .. [1] Chang, S. Grace, Bin Yu, and Martin Vetterli. "Adaptive wavelet
@@ -390,24 +401,20 @@ def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft',
            DOI: 10.1093/biomet/81.3.425
 
     """
-    if not isinstance(wavelet, pywt.Wavelet):
-        wavelet = pywt.Wavelet(wavelet)
+    wavelet = pywt.Wavelet(wavelet)
 
-    # determine the number of wavelet decomposition levels
+    # Determine the number of wavelet decomposition levels
     if wavelet_levels is None:
-        # determine the maximum number of possible levels for img
+        # Determine the maximum number of possible levels for img
         dlen = wavelet.dec_len
         wavelet_levels = np.min(
             [pywt.dwt_max_level(s, dlen) for s in img.shape])
 
-        # BayesShrink variance estimation doesn't work well on levels with
-        # extremely small coefficient arrays, so skip a few of the coarsest
-        # levels.
-        # Note: ref [1] used a fixed wavelet_levels = 4
+        # Skip coarsest wavelet scales (see Notes in docstring).
         wavelet_levels = max(wavelet_levels - 3, 1)
 
     coeffs = pywt.wavedecn(img, wavelet=wavelet, level=wavelet_levels)
-    # detail coefficients at each decomposition level
+    # Detail coefficients at each decomposition level
     dcoeffs = coeffs[1:]
 
     if sigma is None:
@@ -422,13 +429,13 @@ def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft',
                      for level in dcoeffs]
 
     if np.isscalar(threshold):
-        # a single threshold for all coefficient arrays
+        # A single threshold for all coefficient arrays
         denoised_detail = [{key: pywt.threshold(level[key],
                                                 value=threshold,
                                                 mode=mode) for key in level}
                            for level in dcoeffs]
     else:
-        # dict of unique threshold coefficients for each detail coeff. array
+        # Dict of unique threshold coefficients for each detail coeff. array
         denoised_detail = [{key: pywt.threshold(level[key],
                                                 value=thresh[key],
                                                 mode=mode) for key in level}

From 04267966156fcdbf0e7dd7625b1be61c2b49fdbe Mon Sep 17 00:00:00 2001
From: "Gregory R. Lee" <gregory.lee@cchmc.org>
Date: Mon, 15 Aug 2016 13:00:31 -0400
Subject: [PATCH 08/16] TST: test _wavelet_threshold with a single scalar value
 for all levels

---
 skimage/restoration/tests/test_denoise.py | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py
index 6135e6828a9..aad923cd54c 100644
--- a/skimage/restoration/tests/test_denoise.py
+++ b/skimage/restoration/tests/test_denoise.py
@@ -339,6 +339,21 @@ def test_wavelet_denoising():
         assert np.sum(res1**2) <= np.sum(res2**2)
 
 
+def test_wavelet_threshold():
+    rstate = np.random.RandomState(1234)
+
+    img = astro_gray
+    sigma = 0.1
+    noisy = img + sigma * rstate.randn(*(img.shape))
+    noisy = np.clip(noisy, 0, 1)
+
+    # employ a single, uniform threshold instead of BayesShrink sigmas
+    denoised = _wavelet_threshold(noisy, wavelet='db1', threshold=sigma)
+    psnr_noisy = compare_psnr(img, noisy)
+    psnr_denoised = compare_psnr(img, denoised)
+    assert psnr_denoised > psnr_noisy
+
+
 def test_wavelet_denoising_nd():
     rstate = np.random.RandomState(1234)
     for ndim in range(1, 5):

From 915ae055e5ebf6eacaa86daf0c896898dadedcbc Mon Sep 17 00:00:00 2001
From: "Gregory R. Lee" <gregory.lee@cchmc.org>
Date: Mon, 15 Aug 2016 13:42:58 -0400
Subject: [PATCH 09/16] TST: fix call to _wavelet_threshold

---
 skimage/restoration/tests/test_denoise.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py
index aad923cd54c..f76362cd8d3 100644
--- a/skimage/restoration/tests/test_denoise.py
+++ b/skimage/restoration/tests/test_denoise.py
@@ -5,6 +5,7 @@
 from skimage import restoration, data, color, img_as_float, measure
 from skimage._shared._warnings import expected_warnings
 from skimage.measure import compare_psnr
+from skimage.restoration._denoise import _wavelet_threshold
 
 np.random.seed(1234)
 

From b7a70779153e1ce706086b0fc6e524a199660ca9 Mon Sep 17 00:00:00 2001
From: "Gregory R. Lee" <gregory.lee@cchmc.org>
Date: Tue, 6 Sep 2016 17:48:17 -0400
Subject: [PATCH 10/16] ENH: make wavelet_levels an optional parameter to
 denoise_wavelet

---
 skimage/restoration/_denoise.py | 15 +++++++++------
 1 file changed, 9 insertions(+), 6 deletions(-)

diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index 2f9f9fbcc08..e79cd9c3114 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -346,7 +346,7 @@ def _bayes_thresh(details, var):
 
 def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft',
                        wavelet_levels=None):
-    """Performs wavelet denoising.
+    """Perform wavelet denoising.
 
     Parameters
     ----------
@@ -445,8 +445,8 @@ def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft',
 
 
 def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft',
-                    multichannel=False):
-    """Performs wavelet denoising on an image.
+                    multichannel=False, wavelet_levels=None):
+    """Perform wavelet denoising on an image.
 
     Parameters
     ----------
@@ -469,6 +469,9 @@ def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft',
     multichannel : bool, optional
         Apply wavelet denoising separately for each channel (where channels
         correspond to the final axis of the array).
+    wavelet_levels : int or None, optional
+        The number of wavelet decomposition levels to use.  The default is
+        three less than the maximum number of possible decomposition levels.
 
     Returns
     -------
@@ -508,17 +511,17 @@ def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft',
     >>> denoised_img = denoise_wavelet(img, sigma=0.1)
 
     """
-
     img = img_as_float(img)
 
     if multichannel:
         out = np.empty_like(img)
         for c in range(img.shape[-1]):
             out[..., c] = _wavelet_threshold(img[..., c], wavelet=wavelet,
-                                             mode=mode, sigma=sigma)
+                                             mode=mode, sigma=sigma,
+                                             wavelet_levels=wavelet_levels)
     else:
         out = _wavelet_threshold(img, wavelet=wavelet, mode=mode,
-                                 sigma=sigma)
+                                 sigma=sigma, wavelet_levels=wavelet_levels)
 
     clip_range = (-1, 1) if img.min() < 0 else (0, 1)
     return np.clip(out, *clip_range)

From 34d3b59a802e4d1c9e571d1ae6658c08301fbd7f Mon Sep 17 00:00:00 2001
From: "Gregory R. Lee" <gregory.lee@cchmc.org>
Date: Tue, 6 Sep 2016 17:50:35 -0400
Subject: [PATCH 11/16] STY: fix typo in comments

---
 skimage/restoration/_denoise.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index e79cd9c3114..635ae2fbefa 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -418,7 +418,7 @@ def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft',
     dcoeffs = coeffs[1:]
 
     if sigma is None:
-        # Estimates via the noise via method in [2]
+        # Estimate the noise via the method in [2]_
         detail_coeffs = dcoeffs[-1]['d' * img.ndim]
         sigma = np.median(np.abs(detail_coeffs)) / 0.67448975019608171
 

From 3a8d302eedc85d9f42a9ae86701089b785e0e618 Mon Sep 17 00:00:00 2001
From: "Gregory R. Lee" <gregory.lee@cchmc.org>
Date: Tue, 6 Sep 2016 18:13:40 -0400
Subject: [PATCH 12/16] TST: add tests for custom wavelet decomposition levels

---
 skimage/restoration/tests/test_denoise.py | 34 +++++++++++++++++++++++
 1 file changed, 34 insertions(+)

diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py
index f76362cd8d3..a3462595cd0 100644
--- a/skimage/restoration/tests/test_denoise.py
+++ b/skimage/restoration/tests/test_denoise.py
@@ -7,6 +7,8 @@
 from skimage.measure import compare_psnr
 from skimage.restoration._denoise import _wavelet_threshold
 
+import pywt
+
 np.random.seed(1234)
 
 
@@ -373,6 +375,38 @@ def test_wavelet_denoising_nd():
         assert psnr_denoised > psnr_noisy
 
 
+def test_wavelet_denoising_levels():
+    rstate = np.random.RandomState(1234)
+    ndim = 2
+    N = 256
+    wavelet = 'db1'
+    # Generate a very simple test image
+    img = 0.2*np.ones((N, )*ndim)
+    img[[slice(5, 13), ] * ndim] = 0.8
+
+    sigma = 0.1
+    noisy = img + sigma * rstate.randn(*(img.shape))
+    noisy = np.clip(noisy, 0, 1)
+
+    denoised = restoration.denoise_wavelet(noisy, wavelet=wavelet)
+    denoised_1 = restoration.denoise_wavelet(noisy, wavelet=wavelet,
+                                             wavelet_levels=1)
+    psnr_noisy = compare_psnr(img, noisy)
+    psnr_denoised = compare_psnr(img, denoised)
+    psnr_denoised_1 = compare_psnr(img, denoised_1)
+
+    # multi-level case should outperform single level case
+    assert psnr_denoised > psnr_denoised_1 > psnr_noisy
+
+    # invalid number of wavelet levels results in a ValueError
+    max_level = pywt.dwt_max_level(np.min(img.shape),
+                                   pywt.Wavelet(wavelet).dec_len)
+    assert_raises(ValueError, restoration.denoise_wavelet, noisy,
+                  wavelet=wavelet, wavelet_levels=max_level+1)
+    assert_raises(ValueError, restoration.denoise_wavelet, noisy,
+                  wavelet=wavelet, wavelet_levels=-1)
+
+
 def test_estimate_sigma_gray():
     rstate = np.random.RandomState(1234)
     # astronaut image

From 251fe2f4432bfb6bbb83ceca944dde83a7ee2949 Mon Sep 17 00:00:00 2001
From: "Gregory R. Lee" <gregory.lee@cchmc.org>
Date: Wed, 7 Sep 2016 14:35:43 -0400
Subject: [PATCH 13/16] STY: move multichannel argument to last in
 denoise_wavelet

---
 skimage/restoration/_denoise.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index 635ae2fbefa..b68f8dac77a 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -445,7 +445,7 @@ def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft',
 
 
 def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft',
-                    multichannel=False, wavelet_levels=None):
+                    wavelet_levels=None, multichannel=False):
     """Perform wavelet denoising on an image.
 
     Parameters

From deb96a2201036003a83fc8239d542475fe3d9253 Mon Sep 17 00:00:00 2001
From: "Gregory R. Lee" <gregory.lee@cchmc.org>
Date: Wed, 7 Sep 2016 14:37:30 -0400
Subject: [PATCH 14/16] MAINT: reuse code from _sigma_est_dwt within
 _wavelet_threshold

---
 skimage/restoration/_denoise.py | 74 ++++++++++++++++-----------------
 1 file changed, 36 insertions(+), 38 deletions(-)

diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index b68f8dac77a..83b1614b703 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -344,6 +344,41 @@ def _bayes_thresh(details, var):
     return thresh
 
 
+def _sigma_est_dwt(detail_coeffs, distribution='Gaussian'):
+    """Calculate the robust median estimator of the noise standard deviation.
+
+    Parameters
+    ----------
+    detail_coeffs : ndarray
+        The detail coefficients corresponding to the discrete wavelet
+        transform of an image.
+    distribution : str
+        The underlying noise distribution.
+
+    Returns
+    -------
+    sigma : float
+        The estimated noise standard deviation (see section 4.2 of [1]_).
+
+    References
+    ----------
+    .. [1] D. L. Donoho and I. M. Johnstone. "Ideal spatial adaptation
+       by wavelet shrinkage." Biometrika 81.3 (1994): 425-455.
+       DOI:10.1093/biomet/81.3.425
+    """
+    # Consider regions with detail coefficients exactly zero to be masked out
+    detail_coeffs = detail_coeffs[np.nonzero(detail_coeffs)]
+
+    if distribution.lower() == 'gaussian':
+        # 75th quantile of the underlying, symmetric noise distribution
+        denom = scipy.stats.norm.ppf(0.75)
+        sigma = np.median(np.abs(detail_coeffs)) / denom
+    else:
+        raise ValueError("Only Gaussian noise estimation is currently "
+                         "supported")
+    return sigma
+
+
 def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft',
                        wavelet_levels=None):
     """Perform wavelet denoising.
@@ -420,7 +455,7 @@ def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft',
     if sigma is None:
         # Estimate the noise via the method in [2]_
         detail_coeffs = dcoeffs[-1]['d' * img.ndim]
-        sigma = np.median(np.abs(detail_coeffs)) / 0.67448975019608171
+        sigma = _sigma_est_dwt(detail_coeffs, distribution='Gaussian')
 
     if threshold is None:
         # The BayesShrink thresholds from [1]_ in docstring
@@ -527,43 +562,6 @@ def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft',
     return np.clip(out, *clip_range)
 
 
-def _sigma_est_dwt(detail_coeffs, distribution='Gaussian'):
-    """
-    Calculation of the robust median estimator of the noise standard
-    deviation.
-
-    Parameters
-    ----------
-    detail_coeffs : ndarray
-        The detail coefficients corresponding to the discrete wavelet
-        transform of an image.
-    distribution : str
-        The underlying noise distribution.
-
-    Returns
-    -------
-    sigma : float
-        The estimated noise standard deviation (see section 4.2 of [1]_).
-
-    References
-    ----------
-    .. [1] D. L. Donoho and I. M. Johnstone. "Ideal spatial adaptation
-       by wavelet shrinkage." Biometrika 81.3 (1994): 425-455.
-       DOI:10.1093/biomet/81.3.425
-    """
-    # consider regions with detail coefficients exactly zero to be masked out
-    detail_coeffs = detail_coeffs[np.nonzero(detail_coeffs)]
-
-    if distribution.lower() == 'gaussian':
-        # 75th quantile of the underlying, symmetric noise distribution:
-        denom = scipy.stats.norm.ppf(0.75)
-        sigma = np.median(np.abs(detail_coeffs)) / denom
-    else:
-        raise ValueError("Only Gaussian noise estimation is currently "
-                         "supported")
-    return sigma
-
-
 def estimate_sigma(im, multichannel=False, average_sigmas=False):
     """
     Robust wavelet-based estimator of the (Gaussian) noise standard deviation.

From 1e38a55f0d54e9c67abd2c5c6e291f882361a888 Mon Sep 17 00:00:00 2001
From: "Gregory R. Lee" <gregory.lee@cchmc.org>
Date: Wed, 7 Sep 2016 22:04:58 -0400
Subject: [PATCH 15/16] DOC: fix order of parameters in docstring

---
 skimage/restoration/_denoise.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index 83b1614b703..ff412441485 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -501,12 +501,12 @@ def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft',
         An optional argument to choose the type of denoising performed. It
         noted that choosing soft thresholding given additive noise finds the
         best approximation of the original image.
-    multichannel : bool, optional
-        Apply wavelet denoising separately for each channel (where channels
-        correspond to the final axis of the array).
     wavelet_levels : int or None, optional
         The number of wavelet decomposition levels to use.  The default is
         three less than the maximum number of possible decomposition levels.
+    multichannel : bool, optional
+        Apply wavelet denoising separately for each channel (where channels
+        correspond to the final axis of the array).
 
     Returns
     -------

From c4d80582b65366f5c93a71650c1b7a337d43b487 Mon Sep 17 00:00:00 2001
From: "Gregory R. Lee" <gregory.lee@cchmc.org>
Date: Wed, 7 Sep 2016 22:07:07 -0400
Subject: [PATCH 16/16] make multichannel the last argument of estimate_sigma

---
 skimage/restoration/_denoise.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index ff412441485..3c2d46438a6 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -562,7 +562,7 @@ def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft',
     return np.clip(out, *clip_range)
 
 
-def estimate_sigma(im, multichannel=False, average_sigmas=False):
+def estimate_sigma(im, average_sigmas=False, multichannel=False):
     """
     Robust wavelet-based estimator of the (Gaussian) noise standard deviation.
 
@@ -570,11 +570,11 @@ def estimate_sigma(im, multichannel=False, average_sigmas=False):
     ----------
     im : ndarray
         Image for which to estimate the noise standard deviation.
-    multichannel : bool
-        Estimate sigma separately for each channel.
     average_sigmas : bool, optional
         If true, average the channel estimates of `sigma`.  Otherwise return
         a list of sigmas corresponding to each channel.
+    multichannel : bool
+        Estimate sigma separately for each channel.
 
     Returns
     -------