From c83fe7d178a959e87d388f663b29865266784c79 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Wed, 10 Aug 2016 18:32:07 -0400 Subject: [PATCH 1/8] ENH: add nd support to denoise_wavelet --- skimage/restoration/_denoise.py | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py index a853a63c..795d3602 100644 --- a/skimage/restoration/_denoise.py +++ b/skimage/restoration/_denoise.py @@ -394,7 +394,7 @@ def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft'): return pywt.waverecn(denoised_coeffs, wavelet) -def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft'): +def denoise_wavelet(img, multichannel, sigma=None, wavelet='db1', mode='soft'): """Performs wavelet denoising on an image. Parameters @@ -457,16 +457,13 @@ def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft'): img = img_as_float(img) - if img.ndim not in {2, 3}: - raise ValueError('denoise_wavelet only supports 2D and 3D images') - - if img.ndim == 2: + if multichannel: + out = np.stack([_wavelet_threshold(img[..., c], wavelet=wavelet, + mode=mode, sigma=sigma) + for c in range(img.ndim)], axis=-1) + else: out = _wavelet_threshold(img, wavelet=wavelet, mode=mode, sigma=sigma) - else: - out = np.dstack([_wavelet_threshold(img[..., c], wavelet=wavelet, - mode=mode, sigma=sigma) - for c in range(img.ndim)]) clip_range = (-1, 1) if img.min() < 0 else (0, 1) return np.clip(out, *clip_range) From d4168b665ee8735178f1be989bbda4cb314fe1fe Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Wed, 10 Aug 2016 18:33:02 -0400 Subject: [PATCH 2/8] MAINT: more memory efficient multichannel implementation also make multichannel default to False to match other denoising routines --- skimage/restoration/_denoise.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py index 795d3602..4e724447 100644 --- a/skimage/restoration/_denoise.py +++ b/skimage/restoration/_denoise.py @@ -394,7 +394,8 @@ def _wavelet_threshold(img, wavelet, threshold=None, sigma=None, mode='soft'): return pywt.waverecn(denoised_coeffs, wavelet) -def denoise_wavelet(img, multichannel, sigma=None, wavelet='db1', mode='soft'): +def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft', + multichannel=False): """Performs wavelet denoising on an image. Parameters @@ -458,9 +459,10 @@ def denoise_wavelet(img, multichannel, sigma=None, wavelet='db1', mode='soft'): img = img_as_float(img) if multichannel: - out = np.stack([_wavelet_threshold(img[..., c], wavelet=wavelet, - mode=mode, sigma=sigma) - for c in range(img.ndim)], axis=-1) + out = np.empty_like(img) + for c in range(img.ndim): + out[..., c] = _wavelet_threshold(img[..., c], wavelet=wavelet, + mode=mode, sigma=sigma) else: out = _wavelet_threshold(img, wavelet=wavelet, mode=mode, sigma=sigma) From efb829f876a5fdff62445b7cd7f79de280734f3b Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Wed, 10 Aug 2016 18:45:39 -0400 Subject: [PATCH 3/8] TST: update existing tests for wavelet_denoise --- skimage/restoration/tests/test_denoise.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py index 383d0f96..7ee1acd1 100644 --- a/skimage/restoration/tests/test_denoise.py +++ b/skimage/restoration/tests/test_denoise.py @@ -310,16 +310,20 @@ def test_no_denoising_for_small_h(): def test_wavelet_denoising(): - for img in [astro_gray, astro]: + for img, multichannel in [(astro_gray, False), (astro, True)]: noisy = img.copy() + 0.1 * np.random.randn(*(img.shape)) noisy = np.clip(noisy, 0, 1) # less energy in signal - denoised = restoration.denoise_wavelet(noisy, sigma=0.3) + denoised = restoration.denoise_wavelet(noisy, sigma=0.3, + multichannel=multichannel) assert denoised.sum()**2 <= img.sum()**2 # test changing noise_std (higher threshold, so less energy in signal) - assert (restoration.denoise_wavelet(noisy, sigma=0.2).sum()**2 <= - restoration.denoise_wavelet(noisy, sigma=0.1).sum()**2) + res1 = restoration.denoise_wavelet(noisy, sigma=0.2, + multichannel=multichannel) + res2 = restoration.denoise_wavelet(noisy, sigma=0.1, + multichannel=multichannel) + assert (res1.sum()**2 <= res2.sum()**2) if __name__ == "__main__": From 6f3520b83eee7640ea466a8c090106a930f84931 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Thu, 11 Aug 2016 10:51:58 -0400 Subject: [PATCH 4/8] DOC: add missing multichannel info to the docstring --- skimage/restoration/_denoise.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/skimage/restoration/_denoise.py b/skimage/restoration/_denoise.py index 4e724447..377cefc3 100644 --- a/skimage/restoration/_denoise.py +++ b/skimage/restoration/_denoise.py @@ -400,7 +400,7 @@ def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft', Parameters ---------- - img : ndarray (2D/3D) of ints, uints or floats + img : ndarray ([M[, N[, ...P]][, C]) 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. @@ -416,6 +416,9 @@ 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). Returns ------- From 72577e9764fe1cf1920f2f394d0de64c7be71fb7 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Thu, 11 Aug 2016 10:52:37 -0400 Subject: [PATCH 5/8] FIX: correct channel range in the multichannel case --- 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 377cefc3..7f6c6b7b 100644 --- a/skimage/restoration/_denoise.py +++ b/skimage/restoration/_denoise.py @@ -463,7 +463,7 @@ def denoise_wavelet(img, sigma=None, wavelet='db1', mode='soft', if multichannel: out = np.empty_like(img) - for c in range(img.ndim): + for c in range(img.shape[-1]): out[..., c] = _wavelet_threshold(img[..., c], wavelet=wavelet, mode=mode, sigma=sigma) else: From 87936c0b096c57c70ed4b2e897465d218479827f Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Fri, 12 Aug 2016 08:22:18 -0400 Subject: [PATCH 6/8] TST: compare PSNR in the wavelet denoising tests --- skimage/restoration/tests/test_denoise.py | 29 ++++++++++++++++------- 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py index 7ee1acd1..497a472d 100644 --- a/skimage/restoration/tests/test_denoise.py +++ b/skimage/restoration/tests/test_denoise.py @@ -3,7 +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 +from skimage.measure import compare_psnr np.random.seed(1234) @@ -311,17 +311,28 @@ def test_no_denoising_for_small_h(): def test_wavelet_denoising(): for img, multichannel in [(astro_gray, False), (astro, True)]: - noisy = img.copy() + 0.1 * np.random.randn(*(img.shape)) + sigma = 0.1 + noisy = img.copy() + sigma * np.random.randn(*(img.shape)) noisy = np.clip(noisy, 0, 1) - # less energy in signal - denoised = restoration.denoise_wavelet(noisy, sigma=0.3, - multichannel=multichannel) - assert denoised.sum()**2 <= img.sum()**2 - # test changing noise_std (higher threshold, so less energy in signal) - res1 = restoration.denoise_wavelet(noisy, sigma=0.2, + # Verify that SNR is improved when true sigma is used + denoised = restoration.denoise_wavelet(noisy, sigma=sigma, + multichannel=multichannel) + psnr_noisy = compare_psnr(img, noisy) + psnr_denoised = compare_psnr(img, denoised) + assert psnr_denoised > psnr_noisy + + # Verify that SNR is improved with internally estimated sigma + denoised = restoration.denoise_wavelet(noisy, + multichannel=multichannel) + psnr_noisy = compare_psnr(img, noisy) + psnr_denoised = compare_psnr(img, denoised) + assert psnr_denoised > psnr_noisy + + # Test changing noise_std (higher threshold, so less energy in signal) + res1 = restoration.denoise_wavelet(noisy, sigma=2*sigma, multichannel=multichannel) - res2 = restoration.denoise_wavelet(noisy, sigma=0.1, + res2 = restoration.denoise_wavelet(noisy, sigma=sigma, multichannel=multichannel) assert (res1.sum()**2 <= res2.sum()**2) From 48bfe20ed741a50f4ee32d21ab41c99c8787edc8 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Fri, 12 Aug 2016 08:34:43 -0400 Subject: [PATCH 7/8] TST: add basic wavelet denoising tests for 1D through 4D inputs --- skimage/restoration/tests/test_denoise.py | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/skimage/restoration/tests/test_denoise.py b/skimage/restoration/tests/test_denoise.py index 497a472d..acbd44b3 100644 --- a/skimage/restoration/tests/test_denoise.py +++ b/skimage/restoration/tests/test_denoise.py @@ -312,7 +312,7 @@ def test_no_denoising_for_small_h(): def test_wavelet_denoising(): for img, multichannel in [(astro_gray, False), (astro, True)]: sigma = 0.1 - noisy = img.copy() + sigma * np.random.randn(*(img.shape)) + noisy = img + sigma * np.random.randn(*(img.shape)) noisy = np.clip(noisy, 0, 1) # Verify that SNR is improved when true sigma is used @@ -337,5 +337,22 @@ def test_wavelet_denoising(): assert (res1.sum()**2 <= res2.sum()**2) +def test_wavelet_denoising_nd(): + 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 = np.clip(noisy, 0, 1) + + # Verify that SNR is improved with internally estimated sigma + denoised = restoration.denoise_wavelet(noisy) + psnr_noisy = compare_psnr(img, noisy) + psnr_denoised = compare_psnr(img, denoised) + assert psnr_denoised > psnr_noisy + + if __name__ == "__main__": run_module_suite() From 721dd37eea84e1b18e401a4c72ff1dcf828d2282 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Fri, 12 Aug 2016 08:38:50 -0400 Subject: [PATCH 8/8] DOC: update plot_denoise example to use the multichannel argument in denoise_wavelet --- doc/examples/filters/plot_denoise.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/doc/examples/filters/plot_denoise.py b/doc/examples/filters/plot_denoise.py index 9978b818..810aa7e3 100644 --- a/doc/examples/filters/plot_denoise.py +++ b/doc/examples/filters/plot_denoise.py @@ -52,7 +52,8 @@ 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].imshow(denoise_wavelet(noisy, sigma=0.4*astro.std(), + multichannel=True)) ax[0, 3].axis('off') ax[0, 3].set_title('Wavelet') @@ -62,7 +63,8 @@ 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('(more) Bilateral') -ax[1, 3].imshow(denoise_wavelet(noisy, sigma=0.6*astro.std())) +ax[1, 3].imshow(denoise_wavelet(noisy, sigma=0.6*astro.std(), + multichannel=True)) ax[1, 3].axis('off') ax[1, 3].set_title('(more) Wavelet') ax[1, 0].imshow(astro)