From a1cc31e47fc1560a688dc53d7b7aa2d19c16feac Mon Sep 17 00:00:00 2001
From: Scott Sievert <me@scottsievert.com>
Date: Sat, 16 Jul 2016 11:03:43 -0500
Subject: [PATCH 01/12] ENH: implements wavelet denoising

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

diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index 2d9af004..accaeda5 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -332,3 +332,126 @@ def denoise_tv_chambolle(im, weight=0.1, eps=2.e-4, n_iter_max=200,
     else:
         out = _denoise_tv_chambolle_nd(im, weight, eps, n_iter_max)
     return out
+
+
+def _wavelet_threshold(im, wavelet, threshold=None, sigma=None, mode='soft'):
+    """Performs wavelet denoising.
+    Parameters
+    ----------
+    im : ndarray (2d or 3d) of ints, uints or floats
+        Input data to be denoised. `im` can be of any numeric type,
+        but it is cast into an ndarray of floats for the computation
+        of the denoised image.
+    wavelet : string
+        The type of wavelet to perform. Can be any of the options
+        [pywt.wavelist]_ outputs. For example, this may be any of ``{db1, db2,
+        db3, db4, haar}``.
+    sigma : float, optional
+        The standard deviation of the noise. The noise is estimated when sigma is None (the default).
+    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.
+    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
+        best approximation of the original image.
+
+    Returns
+    -------
+    out : ndarray
+        Denoised image.
+
+    References
+    ----------
+    .. [1] Chang, S. Grace, Bin Yu, and Martin Vetterli. "Adaptive wavelet
+           thresholding for image denoising and compression." Image Processing,
+           IEEE Transactions on 9.9 (2000): 1532-1546.
+    """
+    import pywt
+    coeffs = pywt.wavedecn(im, wavelet=wavelet)
+    detail_coeffs = coeffs[-1]['d' * im.ndim]
+
+    if sigma is None:
+        # Estimate the noise std.dev as discussed in PR #1837
+        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(im.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_root = pywt.threshold(coeffs[0], value=threshold, mode=mode)
+    return pywt.waverecn([denoised_root, *denoised_detail], wavelet)
+
+
+def denoise_wavelet(im, sigma=None, wavelet='db1', mode='soft'):
+    """Performs wavelet denoising on an image.
+    Parameters
+    ----------
+    im : ndarray (greater than 2d) of ints, uints or floats
+        Input data to be denoised. `im` can be of any numeric type,
+        but it is cast into an ndarray of floats for the computation
+        of the denoised image.
+    sigma : float, optional
+        The noise standard deviation used when computing the threshold
+        adaptively as described in [1]_.
+    wavelet : string, optional
+        The type of wavelet to perform and can be any of the options
+        [pywt.wavelist]_ outputs. The default is `'db1'`. For example,
+        ``wavelet`` can be any of ``{'db2', 'haar', 'sym9'}`` and many more.
+    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
+        best approximation of the original image.
+    Returns
+    -------
+    out : ndarray
+        Denoised image.
+    Notes
+    -----
+    As with the Fourier transform, there is an analogue to frequency in the
+    wavelet domain. Correspondingly, many pixel values of an image are 0 after
+    taking the wavelet transform.
+    By wavelet denoising, we are enforcing that many of the wavelet coefficients
+    are 0 while keeping the error small. When we use soft thresholding, our
+    estimate is
+    .. math:: \widehat{x} = \arg \min_x ||z - x||_2^2 + \lambda ||x||_1
+    where :math:`z` is the input image wavelet coefficients and :math:`\lambda`
+    is the threshold.
+    This function performs wavelet denoising on each color plane separately. The
+    output is clipped between 0 and 1.
+    References
+    ----------
+    .. [1] Chang, S. Grace, Bin Yu, and Martin Vetterli. "Adaptive wavelet
+           thresholding for image denoising and compression." Image Processing,
+           IEEE Transactions on 9.9 (2000): 1532-1546.
+    .. [pywt.wavelist] http://pywavelets.readthedocs.org/en/latest/ref/wavelets.html#wavelet-wavelist
+
+    Examples
+    --------
+    >>> from skimage import color, data
+    >>> img = data.astronaut() * 1.0 / 255
+    >>> img = color.rgb2gray(img)
+    >>> img += 0.5 * img.std() * np.random.randn(*img.shape)
+    >>> img = np.clip(img, 0, 1)
+    >>> denoised_img = denoise_wavelet(img)
+    >>> assert denoised_img.min() >= 0.0
+    >>> assert denoised_img.max() <= 1.0
+    """
+
+    if not im.dtype.kind == 'f':
+        im = img_as_float(im)
+
+    if im.ndim == 2:
+        out = _wavelet_threshold(im, wavelet=wavelet, mode=mode,
+                                 sigma=sigma)
+
+    else:
+        out = np.dstack([_wavelet_threshold(im[..., c], wavelet=wavelet,
+                                            mode=mode, sigma=sigma)
+                         for c in range(im.ndim)])
+
+    # ensure valid image in 0, 1 is returned
+    return np.clip(out, 0, 1)

From fb8f09263c44fd254634a73730211dd2e19ee22a Mon Sep 17 00:00:00 2001
From: Scott Sievert <me@scottsievert.com>
Date: Sat, 16 Jul 2016 11:03:58 -0500
Subject: [PATCH 02/12] TST: tests wavelet denoising

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

diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py
index 12981d07..8a0fc7de 100644
--- a/skimage/restoration/tests/test_denoise.py
+++ b/skimage/restoration/tests/test_denoise.py
@@ -308,5 +308,21 @@ def test_no_denoising_for_small_h():
     assert np.allclose(denoised, img)
 
 
+def test_wavelet_denoising():
+    img = astro_gray.copy() + 0.1 * np.random.randn(*(astro_gray.shape))
+    img = np.clip(img, 0, 1)
+
+    # less energy in signal
+    assert restoration.denoise_wavelet(img).sum() <= img.sum()
+
+    # test changing noise_std (higher threshold, so less energy in signal)
+    assert (restoration.denoise_wavelet(img, noise_stdev=0.2).sum() <=
+            restoration.denoise_wavelet(img, noise_stdev=0.1).sum())
+
+    # This works for this particular image (probably more wavelet theory here)
+    assert (np.sum(restoration.denoise_wavelet(img, wavelet='db2')) <
+            np.sum(restoration.denoise_wavelet(img, wavelet='db1')))
+
+
 if __name__ == "__main__":
     run_module_suite()

From dfdd0fbf5a76afa7cd8bebb71edc5e3cb90d8234 Mon Sep 17 00:00:00 2001
From: Scott Sievert <me@scottsievert.com>
Date: Sat, 16 Jul 2016 11:04:18 -0500
Subject: [PATCH 03/12] DOC: shows example of wavelet denoising

---
 doc/examples/filters/plot_denoise.py | 24 +++++++++++++++---------
 1 file changed, 15 insertions(+), 9 deletions(-)

diff --git a/doc/examples/filters/plot_denoise.py b/doc/examples/filters/plot_denoise.py
index 7896e656..120b6d11 100644
--- a/doc/examples/filters/plot_denoise.py
+++ b/doc/examples/filters/plot_denoise.py
@@ -29,7 +29,7 @@ import numpy as np
 import matplotlib.pyplot as plt
 
 from skimage import data, img_as_float
-from skimage.restoration import denoise_tv_chambolle, denoise_bilateral
+from skimage.restoration import denoise_tv_chambolle, denoise_bilateral, denoise_wavelet
 
 
 astro = img_as_float(data.astronaut())
@@ -38,7 +38,7 @@ astro = astro[220:300, 220:320]
 noisy = astro + 0.6 * astro.std() * np.random.random(astro.shape)
 noisy = np.clip(noisy, 0, 1)
 
-fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(8, 5), sharex=True,
+fig, ax = plt.subplots(nrows=2, ncols=4, figsize=(8, 5), sharex=True,
                        sharey=True, subplot_kw={'adjustable': 'box-forced'})
 
 plt.gray()
@@ -52,16 +52,22 @@ ax[0, 1].set_title('TV')
 ax[0, 2].imshow(denoise_bilateral(noisy, sigma_color=0.05, sigma_spatial=15))
 ax[0, 2].axis('off')
 ax[0, 2].set_title('Bilateral')
+ax[0, 3].imshow(denoise_wavelet(noisy, sigma=0.4*astro.std()))
+ax[0, 3].axis('off')
+ax[0, 3].set_title('Wavelet')
 
-ax[1, 0].imshow(denoise_tv_chambolle(noisy, weight=0.2, multichannel=True))
-ax[1, 0].axis('off')
-ax[1, 0].set_title('(more) TV')
-ax[1, 1].imshow(denoise_bilateral(noisy, sigma_color=0.1, sigma_spatial=15))
+ax[1, 1].imshow(denoise_tv_chambolle(noisy, weight=0.2, multichannel=True))
 ax[1, 1].axis('off')
-ax[1, 1].set_title('(more) Bilateral')
-ax[1, 2].imshow(astro)
+ax[1, 1].set_title('(more) TV')
+ax[1, 2].imshow(denoise_bilateral(noisy, sigma_color=0.1, sigma_spatial=15))
 ax[1, 2].axis('off')
-ax[1, 2].set_title('original')
+ax[1, 2].set_title('(more) Bilateral')
+ax[1, 3].imshow(denoise_wavelet(noisy, sigma=0.6*astro.std()))
+ax[1, 3].axis('off')
+ax[1, 3].set_title('(more) Wavelet')
+ax[1, 0].imshow(astro)
+ax[1, 0].axis('off')
+ax[1, 0].set_title('original')
 
 fig.tight_layout()
 

From 1c38da0fca2c8b7a461f6f139286996d9484dd4b Mon Sep 17 00:00:00 2001
From: Scott Sievert <me@scottsievert.com>
Date: Sat, 16 Jul 2016 11:07:22 -0500
Subject: [PATCH 04/12] MAINT: adds pywt to optional requirements

---
 optional_requirements.txt | 1 +
 1 file changed, 1 insertion(+)

diff --git a/optional_requirements.txt b/optional_requirements.txt
index 8bc53bb0..f2611edb 100644
--- a/optional_requirements.txt
+++ b/optional_requirements.txt
@@ -4,3 +4,4 @@ SimpleITK; python_version != '3.4' and python_version != '3.5'
 astropy
 tifffile
 imageio
+PyWavelet

From 643fe7784039922e32196ffd6e01acb3ec6c3e8b Mon Sep 17 00:00:00 2001
From: Scott Sievert <me@scottsievert.com>
Date: Sat, 16 Jul 2016 11:51:58 -0500
Subject: [PATCH 05/12] MAINT: first pass features (see below)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

im —> img
adds paper DOI
adds denoise_wavelet __init__
fixes typo in optional_requriments
---
 optional_requirements.txt       |  2 +-
 skimage/restoration/__init__.py |  3 +-
 skimage/restoration/_denoise.py | 57 +++++++++++++++++++--------------
 3 files changed, 36 insertions(+), 26 deletions(-)

diff --git a/optional_requirements.txt b/optional_requirements.txt
index f2611edb..6f923675 100644
--- a/optional_requirements.txt
+++ b/optional_requirements.txt
@@ -4,4 +4,4 @@ SimpleITK; python_version != '3.4' and python_version != '3.5'
 astropy
 tifffile
 imageio
-PyWavelet
+PyWavelets
diff --git a/skimage/restoration/__init__.py b/skimage/restoration/__init__.py
index 424f9939..8e37fb6e 100644
--- a/skimage/restoration/__init__.py
+++ b/skimage/restoration/__init__.py
@@ -21,7 +21,7 @@ References
 from .deconvolution import wiener, unsupervised_wiener, richardson_lucy
 from .unwrap import unwrap_phase
 from ._denoise import denoise_tv_chambolle, denoise_tv_bregman, \
-                      denoise_bilateral
+                      denoise_bilateral, denoise_wavelet
 from .non_local_means import denoise_nl_means
 from .inpaint import inpaint_biharmonic
 from .._shared.utils import copy_func, deprecated
@@ -37,6 +37,7 @@ __all__ = ['wiener',
            'denoise_tv_bregman',
            'denoise_tv_chambolle',
            'denoise_bilateral',
+           'denoise_wavelet',
            'denoise_nl_means',
            'nl_means_denoising',
            'inpaint_biharmonic']
diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index accaeda5..57bc67fd 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -334,11 +334,11 @@ def denoise_tv_chambolle(im, weight=0.1, eps=2.e-4, n_iter_max=200,
     return out
 
 
-def _wavelet_threshold(im, wavelet, threshold=None, sigma=None, mode='soft'):
+def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft'):
     """Performs wavelet denoising.
     Parameters
     ----------
-    im : ndarray (2d or 3d) of ints, uints or floats
+    img : ndarray (2d or 3d) of ints, uints or floats
         Input data to be denoised. `im` can be of any numeric type,
         but it is cast into an ndarray of floats for the computation
         of the denoised image.
@@ -347,11 +347,12 @@ def _wavelet_threshold(im, wavelet, threshold=None, sigma=None, mode='soft'):
         [pywt.wavelist]_ outputs. For example, this may be any of ``{db1, db2,
         db3, db4, haar}``.
     sigma : float, optional
-        The standard deviation of the noise. The noise is estimated when sigma is None (the default).
+        The standard deviation of the noise. The noise is estimated when sigma
+        is None (the default).
     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 SureShrink 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
@@ -367,10 +368,11 @@ def _wavelet_threshold(im, wavelet, threshold=None, sigma=None, mode='soft'):
     .. [1] Chang, S. Grace, Bin Yu, and Martin Vetterli. "Adaptive wavelet
            thresholding for image denoising and compression." Image Processing,
            IEEE Transactions on 9.9 (2000): 1532-1546.
+           DOI: 10.1109/83.862633
     """
     import pywt
-    coeffs = pywt.wavedecn(im, wavelet=wavelet)
-    detail_coeffs = coeffs[-1]['d' * im.ndim]
+    coeffs = pywt.wavedecn(img, wavelet=wavelet)
+    detail_coeffs = coeffs[-1]['d' * img.ndim]
 
     if sigma is None:
         # Estimate the noise std.dev as discussed in PR #1837
@@ -378,7 +380,7 @@ def _wavelet_threshold(im, wavelet, threshold=None, sigma=None, mode='soft'):
 
     if threshold is None:
         # The BayesShrink threshold from [1]_ in docstring
-        threshold = sigma**2 / np.sqrt(max(im.var() - sigma**2, 0))
+        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:]]
@@ -386,11 +388,11 @@ def _wavelet_threshold(im, wavelet, threshold=None, sigma=None, mode='soft'):
     return pywt.waverecn([denoised_root, *denoised_detail], wavelet)
 
 
-def denoise_wavelet(im, sigma=None, wavelet='db1', mode='soft'):
-    """Performs wavelet denoising on an image.
+def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft'):
+    r"""Performs wavelet denoising on an image.
     Parameters
     ----------
-    im : ndarray (greater than 2d) of ints, uints or floats
+    img : ndarray (greater than 2d) of ints, uints or floats
         Input data to be denoised. `im` can be of any numeric type,
         but it is cast into an ndarray of floats for the computation
         of the denoised image.
@@ -405,34 +407,41 @@ def denoise_wavelet(im, 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.
+
     Returns
     -------
     out : ndarray
         Denoised image.
+
     Notes
     -----
     As with the Fourier transform, there is an analogue to frequency in the
     wavelet domain. Correspondingly, many pixel values of an image are 0 after
     taking the wavelet transform.
-    By wavelet denoising, we are enforcing that many of the wavelet coefficients
-    are 0 while keeping the error small. When we use soft thresholding, our
-    estimate is
+
+    By wavelet denoising, we are enforcing that many of the wavelet
+    coefficients are 0 while keeping the error small. When we use soft
+    thresholding, our estimate is
+
     .. math:: \widehat{x} = \arg \min_x ||z - x||_2^2 + \lambda ||x||_1
+
     where :math:`z` is the input image wavelet coefficients and :math:`\lambda`
     is the threshold.
-    This function performs wavelet denoising on each color plane separately. The
-    output is clipped between 0 and 1.
+
+    This function performs wavelet denoising on each color plane separately.
+    The output is clipped between 0 and 1.
+
     References
     ----------
     .. [1] Chang, S. Grace, Bin Yu, and Martin Vetterli. "Adaptive wavelet
            thresholding for image denoising and compression." Image Processing,
            IEEE Transactions on 9.9 (2000): 1532-1546.
-    .. [pywt.wavelist] http://pywavelets.readthedocs.org/en/latest/ref/wavelets.html#wavelet-wavelist
+           DOI: 10.1109/83.862633
 
     Examples
     --------
     >>> from skimage import color, data
-    >>> img = data.astronaut() * 1.0 / 255
+    >>> img = img_as_float(data.astronaut())
     >>> img = color.rgb2gray(img)
     >>> img += 0.5 * img.std() * np.random.randn(*img.shape)
     >>> img = np.clip(img, 0, 1)
@@ -441,17 +450,17 @@ def denoise_wavelet(im, sigma=None, wavelet='db1', mode='soft'):
     >>> assert denoised_img.max() <= 1.0
     """
 
-    if not im.dtype.kind == 'f':
-        im = img_as_float(im)
+    if not img.dtype.kind == 'f':
+        img = img_as_float(img)
 
-    if im.ndim == 2:
-        out = _wavelet_threshold(im, wavelet=wavelet, mode=mode,
+    if img.ndim == 2:
+        out = _wavelet_threshold(img, wavelet=wavelet, mode=mode,
                                  sigma=sigma)
 
     else:
-        out = np.dstack([_wavelet_threshold(im[..., c], wavelet=wavelet,
+        out = np.dstack([_wavelet_threshold(img[..., c], wavelet=wavelet,
                                             mode=mode, sigma=sigma)
-                         for c in range(im.ndim)])
+                         for c in range(img.ndim)])
 
     # ensure valid image in 0, 1 is returned
     return np.clip(out, 0, 1)

From 5300db1e8d1083ded7b7d22882dd46d70eed78bb Mon Sep 17 00:00:00 2001
From: Scott Sievert <me@scottsievert.com>
Date: Sat, 16 Jul 2016 12:09:52 -0500
Subject: [PATCH 06/12] TST: 2d and 3d. DOC: adds refs to doc

---
 skimage/restoration/_denoise.py           | 15 +++++++++++----
 skimage/restoration/tests/test_denoise.py | 22 +++++++++++-----------
 2 files changed, 22 insertions(+), 15 deletions(-)

diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index 57bc67fd..a6bb0047 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -344,11 +344,11 @@ def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft'):
         of the denoised image.
     wavelet : string
         The type of wavelet to perform. Can be any of the options
-        [pywt.wavelist]_ outputs. For example, this may be any of ``{db1, db2,
+        pywt.wavelist outputs. For example, this may be any of ``{db1, db2,
         db3, db4, haar}``.
     sigma : float, optional
         The standard deviation of the noise. The noise is estimated when sigma
-        is None (the default).
+        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
@@ -369,13 +369,16 @@ def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft'):
            thresholding for image denoising and compression." Image Processing,
            IEEE Transactions on 9.9 (2000): 1532-1546.
            DOI: 10.1109/83.862633
+    .. [2] 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
     """
     import pywt
     coeffs = pywt.wavedecn(img, wavelet=wavelet)
     detail_coeffs = coeffs[-1]['d' * img.ndim]
 
     if sigma is None:
-        # Estimate the noise std.dev as discussed in PR #1837
+        # Estimates via the noise via method in [2]
         sigma = np.median(np.abs(detail_coeffs)) / 0.67448975019608171
 
     if threshold is None:
@@ -398,7 +401,8 @@ def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft'):
         of the denoised image.
     sigma : float, optional
         The noise standard deviation used when computing the threshold
-        adaptively as described in [1]_.
+        adaptively as described in [1]_. When None (default), the noise
+        standard deviation is estimated via the method in [2]_.
     wavelet : string, optional
         The type of wavelet to perform and can be any of the options
         [pywt.wavelist]_ outputs. The default is `'db1'`. For example,
@@ -437,6 +441,9 @@ def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft'):
            thresholding for image denoising and compression." Image Processing,
            IEEE Transactions on 9.9 (2000): 1532-1546.
            DOI: 10.1109/83.862633
+    .. [2] 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
 
     Examples
     --------
diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py
index 8a0fc7de..664f89a3 100644
--- a/skimage/restoration/tests/test_denoise.py
+++ b/skimage/restoration/tests/test_denoise.py
@@ -309,19 +309,19 @@ def test_no_denoising_for_small_h():
 
 
 def test_wavelet_denoising():
-    img = astro_gray.copy() + 0.1 * np.random.randn(*(astro_gray.shape))
-    img = np.clip(img, 0, 1)
+    for img in [astro_gray, astro]:
+        img = img.copy() + 0.1 * np.random.randn(*(img.shape))
+        img = np.clip(img, 0, 1)
+        # less energy in signal
+        assert restoration.denoise_wavelet(img).sum() <= img.sum()
 
-    # less energy in signal
-    assert restoration.denoise_wavelet(img).sum() <= img.sum()
+        # test changing noise_std (higher threshold, so less energy in signal)
+        assert (restoration.denoise_wavelet(img, noise_stdev=0.2).sum() <=
+                restoration.denoise_wavelet(img, noise_stdev=0.1).sum())
 
-    # test changing noise_std (higher threshold, so less energy in signal)
-    assert (restoration.denoise_wavelet(img, noise_stdev=0.2).sum() <=
-            restoration.denoise_wavelet(img, noise_stdev=0.1).sum())
-
-    # This works for this particular image (probably more wavelet theory here)
-    assert (np.sum(restoration.denoise_wavelet(img, wavelet='db2')) <
-            np.sum(restoration.denoise_wavelet(img, wavelet='db1')))
+        # This works for this particular image
+        assert (np.sum(restoration.denoise_wavelet(img, wavelet='db2')) <
+                np.sum(restoration.denoise_wavelet(img, wavelet='db1')))
 
 
 if __name__ == "__main__":

From 9616637fe64f78026b20d28f4063cc3d36f0cdb4 Mon Sep 17 00:00:00 2001
From: Scott Sievert <me@scottsievert.com>
Date: Sat, 16 Jul 2016 15:08:55 -0500
Subject: [PATCH 07/12] MAINT: py27 compatibility (no splat in list for
 _wavelet_threshold)

MAINT: addresses blank lines, typos
---
 skimage/restoration/_denoise.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index a6bb0047..33fac776 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -339,7 +339,7 @@ def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft'):
     Parameters
     ----------
     img : ndarray (2d or 3d) of ints, uints or floats
-        Input data to be denoised. `im` can be of any numeric type,
+        Input data to be denoised. `img` can be of any numeric type,
         but it is cast into an ndarray of floats for the computation
         of the denoised image.
     wavelet : string
@@ -388,7 +388,8 @@ 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)
-    return pywt.waverecn([denoised_root, *denoised_detail], wavelet)
+    denoised_coeffs = [denoised_root] + [d for d in denoised_detail]
+    return pywt.waverecn(denoised_coeffs, wavelet)
 
 
 def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft'):
@@ -396,7 +397,7 @@ def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft'):
     Parameters
     ----------
     img : ndarray (greater than 2d) of ints, uints or floats
-        Input data to be denoised. `im` can be of any numeric type,
+        Input data to be denoised. `img` can be of any numeric type,
         but it is cast into an ndarray of floats for the computation
         of the denoised image.
     sigma : float, optional
@@ -463,7 +464,6 @@ def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft'):
     if img.ndim == 2:
         out = _wavelet_threshold(img, wavelet=wavelet, mode=mode,
                                  sigma=sigma)
-
     else:
         out = np.dstack([_wavelet_threshold(img[..., c], wavelet=wavelet,
                                             mode=mode, sigma=sigma)

From c6e5c7095350e931709acd894a54b32c42375822 Mon Sep 17 00:00:00 2001
From: Scott Sievert <me@scottsievert.com>
Date: Sat, 23 Jul 2016 15:09:12 -0500
Subject: [PATCH 08/12] makes PyWavelets required

---
 optional_requirements.txt                 |  1 -
 requirements.txt                          |  1 +
 skimage/restoration/_denoise.py           | 25 ++++++++---------------
 skimage/restoration/tests/test_denoise.py |  6 +++---
 4 files changed, 13 insertions(+), 20 deletions(-)

diff --git a/optional_requirements.txt b/optional_requirements.txt
index 6f923675..8bc53bb0 100644
--- a/optional_requirements.txt
+++ b/optional_requirements.txt
@@ -4,4 +4,3 @@ SimpleITK; python_version != '3.4' and python_version != '3.5'
 astropy
 tifffile
 imageio
-PyWavelets
diff --git a/requirements.txt b/requirements.txt
index 1e6921be..b05402dc 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -5,3 +5,4 @@ six>=1.7.3
 networkx>=1.8
 pillow>=2.1.0
 dask[array]>=0.5.0
+PyWavelets>=0.4.0
diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index 33fac776..94a11608 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -5,6 +5,7 @@ from .. import img_as_float
 from ..restoration._denoise_cy import _denoise_bilateral, _denoise_tv_bregman
 from .._shared.utils import skimage_deprecation, warn
 import warnings
+import pywt
 
 
 def denoise_bilateral(image, win_size=None, sigma_color=None, sigma_spatial=1,
@@ -336,6 +337,7 @@ def denoise_tv_chambolle(im, weight=0.1, eps=2.e-4, n_iter_max=200,
 
 def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft'):
     """Performs wavelet denoising.
+
     Parameters
     ----------
     img : ndarray (2d or 3d) of ints, uints or floats
@@ -372,8 +374,8 @@ def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft'):
     .. [2] 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
+
     """
-    import pywt
     coeffs = pywt.wavedecn(img, wavelet=wavelet)
     detail_coeffs = coeffs[-1]['d' * img.ndim]
 
@@ -393,7 +395,8 @@ def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft'):
 
 
 def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft'):
-    r"""Performs wavelet denoising on an image.
+    """Performs wavelet denoising on an image.
+
     Parameters
     ----------
     img : ndarray (greater than 2d) of ints, uints or floats
@@ -406,7 +409,7 @@ def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft'):
         standard deviation is estimated via the method in [2]_.
     wavelet : string, optional
         The type of wavelet to perform and can be any of the options
-        [pywt.wavelist]_ outputs. The default is `'db1'`. For example,
+        ``pywt.wavelist`` outputs. The default is `'db1'`. For example,
         ``wavelet`` can be any of ``{'db2', 'haar', 'sym9'}`` and many more.
     mode : {'soft', 'hard'}, optional
         An optional argument to choose the type of denoising performed. It
@@ -424,15 +427,6 @@ def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft'):
     wavelet domain. Correspondingly, many pixel values of an image are 0 after
     taking the wavelet transform.
 
-    By wavelet denoising, we are enforcing that many of the wavelet
-    coefficients are 0 while keeping the error small. When we use soft
-    thresholding, our estimate is
-
-    .. math:: \widehat{x} = \arg \min_x ||z - x||_2^2 + \lambda ||x||_1
-
-    where :math:`z` is the input image wavelet coefficients and :math:`\lambda`
-    is the threshold.
-
     This function performs wavelet denoising on each color plane separately.
     The output is clipped between 0 and 1.
 
@@ -451,11 +445,10 @@ def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft'):
     >>> from skimage import color, data
     >>> img = img_as_float(data.astronaut())
     >>> img = color.rgb2gray(img)
-    >>> img += 0.5 * img.std() * np.random.randn(*img.shape)
+    >>> img += 0.1 * np.random.randn(*img.shape)
     >>> img = np.clip(img, 0, 1)
-    >>> denoised_img = denoise_wavelet(img)
-    >>> assert denoised_img.min() >= 0.0
-    >>> assert denoised_img.max() <= 1.0
+    >>> denoised_img = denoise_wavelet(img, sigma=0.1)
+
     """
 
     if not img.dtype.kind == 'f':
diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py
index 664f89a3..206f1923 100644
--- a/skimage/restoration/tests/test_denoise.py
+++ b/skimage/restoration/tests/test_denoise.py
@@ -313,11 +313,11 @@ def test_wavelet_denoising():
         img = img.copy() + 0.1 * np.random.randn(*(img.shape))
         img = np.clip(img, 0, 1)
         # less energy in signal
-        assert restoration.denoise_wavelet(img).sum() <= img.sum()
+        assert restoration.denoise_wavelet(img, sigma=0.3).sum() <= img.sum()
 
         # test changing noise_std (higher threshold, so less energy in signal)
-        assert (restoration.denoise_wavelet(img, noise_stdev=0.2).sum() <=
-                restoration.denoise_wavelet(img, noise_stdev=0.1).sum())
+        assert (restoration.denoise_wavelet(img, sigma=0.2).sum() <=
+                restoration.denoise_wavelet(img, sigma=0.1).sum())
 
         # This works for this particular image
         assert (np.sum(restoration.denoise_wavelet(img, wavelet='db2')) <

From 7cbd2d36a43470baf32ee19c35227573981c216d Mon Sep 17 00:00:00 2001
From: Scott Sievert <me@scottsievert.com>
Date: Wed, 3 Aug 2016 15:03:11 -0500
Subject: [PATCH 09/12] MAINT: uses random_noise in examples

---
 doc/examples/filters/plot_denoise.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/doc/examples/filters/plot_denoise.py b/doc/examples/filters/plot_denoise.py
index 120b6d11..9978b818 100644
--- a/doc/examples/filters/plot_denoise.py
+++ b/doc/examples/filters/plot_denoise.py
@@ -30,13 +30,13 @@ import matplotlib.pyplot as plt
 
 from skimage import data, img_as_float
 from skimage.restoration import denoise_tv_chambolle, denoise_bilateral, denoise_wavelet
+from skimage.util import random_noise
 
 
 astro = img_as_float(data.astronaut())
 astro = astro[220:300, 220:320]
 
-noisy = astro + 0.6 * astro.std() * np.random.random(astro.shape)
-noisy = np.clip(noisy, 0, 1)
+noisy = random_noise(astro, var=(0.6 * astro.std())**2)
 
 fig, ax = plt.subplots(nrows=2, ncols=4, figsize=(8, 5), sharex=True,
                        sharey=True, subplot_kw={'adjustable': 'box-forced'})

From 222a1b3a78a1a562327f8dc5c1e322e0b3abe5a4 Mon Sep 17 00:00:00 2001
From: Scott Sievert <me@scottsievert.com>
Date: Thu, 4 Aug 2016 09:48:35 -0500
Subject: [PATCH 10/12] MAINT: asserts image 2D or 3D, proper clipping

as well as DOC: clarify 2D/3D + rewrite notes for denoise_wavelet
---
 skimage/restoration/_denoise.py | 26 ++++++++++++++++----------
 1 file changed, 16 insertions(+), 10 deletions(-)

diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index 94a11608..4d44ab2f 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -399,7 +399,7 @@ def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft'):
 
     Parameters
     ----------
-    img : ndarray (greater than 2d) of ints, uints or floats
+    img : ndarray (2D/3D) of ints, uints or floats
         Input data to be denoised. `img` can be of any numeric type,
         but it is cast into an ndarray of floats for the computation
         of the denoised image.
@@ -423,12 +423,16 @@ def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft'):
 
     Notes
     -----
-    As with the Fourier transform, there is an analogue to frequency in the
-    wavelet domain. Correspondingly, many pixel values of an image are 0 after
-    taking the wavelet transform.
+    The wavelet domain is a sparse representation of the image, and can be
+    thought of similarly to the frequency domain of the Fourier transform.
+    Sparse representations have most values zero or near-zero and truly random
+    noise is (usually) represented by many small values in the wavelet domain.
+    Setting all values below some threshold to 0 reduces the noise in the
+    image, but larger thresholds also decrease the detail present in the image.
 
-    This function performs wavelet denoising on each color plane separately.
-    The output is clipped between 0 and 1.
+    If the input is 3D, this function performs wavelet denoising on each color
+    plane separately. The output image is clipped between either [-1, 1] and
+    [0, 1] depending on the input image range.
 
     References
     ----------
@@ -451,8 +455,10 @@ def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft'):
 
     """
 
-    if not img.dtype.kind == 'f':
-        img = img_as_float(img)
+    img = img_as_float(img)
+
+    if img.ndims not in {2, 3}:
+        raise ValueError('denoise_wavelet only supports 2D and 3D images')
 
     if img.ndim == 2:
         out = _wavelet_threshold(img, wavelet=wavelet, mode=mode,
@@ -462,5 +468,5 @@ def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft'):
                                             mode=mode, sigma=sigma)
                          for c in range(img.ndim)])
 
-    # ensure valid image in 0, 1 is returned
-    return np.clip(out, 0, 1)
+    clip_range = (-1, 1) if img.min() < 0 else (0, 1)
+    return np.clip(out, *clip_range)

From ce3a3da595686cc1fdc46062dafc3a62bac49a40 Mon Sep 17 00:00:00 2001
From: Scott Sievert <me@scottsievert.com>
Date: Thu, 4 Aug 2016 09:48:58 -0500
Subject: [PATCH 11/12] adds to contributors.txt

---
 CONTRIBUTORS.txt | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index 44214d85..0d90a330 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -237,3 +237,6 @@
 
 - Abdeali Kothari
   Alpha blending to convert from rgba to rgb
+
+- Scott Sievert
+  Wavelet denoising

From 13490a64acc3df90f07719d3faaebdcfe1479731 Mon Sep 17 00:00:00 2001
From: Scott Sievert <me@scottsievert.com>
Date: Thu, 4 Aug 2016 09:57:04 -0500
Subject: [PATCH 12/12] TST MAINT: more clearly represent energy in signal TST:
 removes tests to particular image

---
 skimage/restoration/_denoise.py           |  2 +-
 skimage/restoration/tests/test_denoise.py | 16 +++++++---------
 2 files changed, 8 insertions(+), 10 deletions(-)

diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py
index 4d44ab2f..a853a63c 100644
--- a/skimage/restoration/_denoise.py
+++ b/skimage/restoration/_denoise.py
@@ -457,7 +457,7 @@ def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft'):
 
     img = img_as_float(img)
 
-    if img.ndims not in {2, 3}:
+    if img.ndim not in {2, 3}:
         raise ValueError('denoise_wavelet only supports 2D and 3D images')
 
     if img.ndim == 2:
diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py
index 206f1923..383d0f96 100644
--- a/skimage/restoration/tests/test_denoise.py
+++ b/skimage/restoration/tests/test_denoise.py
@@ -3,6 +3,7 @@ from numpy.testing import run_module_suite, assert_raises, assert_equal
 
 from skimage import restoration, data, color, img_as_float, measure
 from skimage._shared._warnings import expected_warnings
+from skimage.measure import compare_ssim
 
 np.random.seed(1234)
 
@@ -310,18 +311,15 @@ def test_no_denoising_for_small_h():
 
 def test_wavelet_denoising():
     for img in [astro_gray, astro]:
-        img = img.copy() + 0.1 * np.random.randn(*(img.shape))
-        img = np.clip(img, 0, 1)
+        noisy = img.copy() + 0.1 * np.random.randn(*(img.shape))
+        noisy = np.clip(noisy, 0, 1)
         # less energy in signal
-        assert restoration.denoise_wavelet(img, sigma=0.3).sum() <= img.sum()
+        denoised = restoration.denoise_wavelet(noisy, sigma=0.3)
+        assert denoised.sum()**2 <= img.sum()**2
 
         # test changing noise_std (higher threshold, so less energy in signal)
-        assert (restoration.denoise_wavelet(img, sigma=0.2).sum() <=
-                restoration.denoise_wavelet(img, sigma=0.1).sum())
-
-        # This works for this particular image
-        assert (np.sum(restoration.denoise_wavelet(img, wavelet='db2')) <
-                np.sum(restoration.denoise_wavelet(img, wavelet='db1')))
+        assert (restoration.denoise_wavelet(noisy, sigma=0.2).sum()**2 <=
+                restoration.denoise_wavelet(noisy, sigma=0.1).sum()**2)
 
 
 if __name__ == "__main__":