1"""
2====================
3Denoising a picture
4====================
5
6In this example, we denoise a noisy version of a picture using the total
7variation, bilateral, and wavelet denoising filters.
8
9Total variation and bilateral algorithms typically produce "posterized" images
10with flat domains separated by sharp edges. It is possible to change the degree
11of posterization by controlling the tradeoff between denoising and faithfulness
12to the original image.
13
14Total variation filter
15----------------------
16
17The result of this filter is an image that has a minimal total variation norm,
18while being as close to the initial image as possible. The total variation is
19the L1 norm of the gradient of the image.
20
21Bilateral filter
22----------------
23
24A bilateral filter is an edge-preserving and noise reducing filter. It averages
25pixels based on their spatial closeness and radiometric similarity.
26
27Wavelet denoising filter
28------------------------
29
30A wavelet denoising filter relies on the wavelet representation of the image.
31The noise is represented by small values in the wavelet domain which are set to
320.
33
34In color images, wavelet denoising is typically done in the `YCbCr color
35space`_ as denoising in separate color channels may lead to more apparent
36noise.
37
38.. _`YCbCr color space`: https://en.wikipedia.org/wiki/YCbCr
39
40"""
41import matplotlib.pyplot as plt
42
43from skimage.restoration import (denoise_tv_chambolle, denoise_bilateral,
44                                 denoise_wavelet, estimate_sigma)
45from skimage import data, img_as_float
46from skimage.util import random_noise
47
48
49original = img_as_float(data.chelsea()[100:250, 50:300])
50
51sigma = 0.155
52noisy = random_noise(original, var=sigma**2)
53
54fig, ax = plt.subplots(nrows=2, ncols=4, figsize=(8, 5),
55                       sharex=True, sharey=True)
56
57plt.gray()
58
59# Estimate the average noise standard deviation across color channels.
60sigma_est = estimate_sigma(noisy, channel_axis=-1, average_sigmas=True)
61# Due to clipping in random_noise, the estimate will be a bit smaller than the
62# specified sigma.
63print(f'Estimated Gaussian noise standard deviation = {sigma_est}')
64
65ax[0, 0].imshow(noisy)
66ax[0, 0].axis('off')
67ax[0, 0].set_title('Noisy')
68ax[0, 1].imshow(denoise_tv_chambolle(noisy, weight=0.1, channel_axis=-1))
69ax[0, 1].axis('off')
70ax[0, 1].set_title('TV')
71ax[0, 2].imshow(denoise_bilateral(noisy, sigma_color=0.05, sigma_spatial=15,
72                channel_axis=-1))
73ax[0, 2].axis('off')
74ax[0, 2].set_title('Bilateral')
75ax[0, 3].imshow(denoise_wavelet(noisy, channel_axis=-1, rescale_sigma=True))
76ax[0, 3].axis('off')
77ax[0, 3].set_title('Wavelet denoising')
78
79ax[1, 1].imshow(denoise_tv_chambolle(noisy, weight=0.2, channel_axis=-1))
80ax[1, 1].axis('off')
81ax[1, 1].set_title('(more) TV')
82ax[1, 2].imshow(denoise_bilateral(noisy, sigma_color=0.1, sigma_spatial=15,
83                channel_axis=-1))
84ax[1, 2].axis('off')
85ax[1, 2].set_title('(more) Bilateral')
86ax[1, 3].imshow(denoise_wavelet(noisy, channel_axis=-1, convert2ycbcr=True,
87                                rescale_sigma=True))
88ax[1, 3].axis('off')
89ax[1, 3].set_title('Wavelet denoising\nin YCbCr colorspace')
90ax[1, 0].imshow(original)
91ax[1, 0].axis('off')
92ax[1, 0].set_title('Original')
93
94fig.tight_layout()
95
96plt.show()
97