pywt

1
"""Using the FSWT to process anistropic images.
2

3
In this demo, an anisotropic piecewise-constant image is transformed by the
4
standard DWT and the fully-separable DWT. The 'Haar' wavelet gives a sparse
5
representation for such piecewise constant signals (detail coefficients are
6
only non-zero near edges).
7

8
For such anistropic signals, the number of non-zero coefficients will be lower
9
for the fully separable DWT than for the isotropic one.
10

11
This example is inspired by the following publication where it is proven that
12
the FSWT gives a sparser representation than the DWT for this class of
13
anistropic images:
14

15
.. V Velisavljevic, B Beferull-Lozano, M Vetterli and PL Dragotti.
16
   Directionlets: Anisotropic Multidirectional Representation With
17
   Separable Filtering. IEEE Transactions on Image Processing, Vol. 15,
18
   No. 7, July 2006.
19

20
"""
21

22
import numpy as np
23
from matplotlib import pyplot as plt
24

25
import pywt
26

27

28
def mondrian(shape=(256, 256), nx=5, ny=8, seed=4):
29
    """ Piecewise-constant image (reminiscent of Dutch painter Piet Mondrian's
30
    geometrical period).
31
    """
32
    rstate = np.random.RandomState(seed)
33
    min_dx = 0
34
    while(min_dx < 3):
35
        xp = np.sort(np.round(rstate.rand(nx-1)*shape[0]).astype(np.int64))
36
        xp = np.concatenate(((0, ), xp, (shape[0], )))
37
        min_dx = np.min(np.diff(xp))
38
    min_dy = 0
39
    while(min_dy < 3):
40
        yp = np.sort(np.round(rstate.rand(ny-1)*shape[1]).astype(np.int64))
41
        yp = np.concatenate(((0, ), yp, (shape[1], )))
42
        min_dy = np.min(np.diff(yp))
43
    img = np.zeros(shape)
44
    for ix, x in enumerate(xp[:-1]):
45
        for iy, y in enumerate(yp[:-1]):
46
            slices = [slice(x, xp[ix+1]), slice(y, yp[iy+1])]
47
            val = rstate.rand(1)[0]
48
            img[slices] = val
49
    return img
50

51

52
# create an anisotropic piecewise constant image
53
img = mondrian((128, 128))
54

55
# perform DWT
56
coeffs_dwt = pywt.wavedecn(img, wavelet='db1', level=None)
57

58
# convert coefficient dictionary to a single array
59
coeff_array_dwt, _ = pywt.coeffs_to_array(coeffs_dwt)
60

61
# perform fully separable DWT
62
fswavedecn_result = pywt.fswavedecn(img, wavelet='db1')
63

64
nnz_dwt = np.sum(coeff_array_dwt != 0)
65
nnz_fswavedecn = np.sum(fswavedecn_result.coeffs != 0)
66

67
print(f"Number of nonzero wavedecn coefficients = {np.sum(nnz_dwt)}")
68
print(f"Number of nonzero fswavedecn coefficients = {np.sum(nnz_fswavedecn)}")
69

70
img = mondrian()
71
fig, axes = plt.subplots(1, 3)
72
imshow_kwargs = {'cmap': plt.cm.gray, 'interpolation': 'nearest'}
73
axes[0].imshow(img, **imshow_kwargs)
74
axes[0].set_title('Anisotropic Image')
75
axes[1].imshow(coeff_array_dwt != 0, **imshow_kwargs)
76
axes[1].set_title(f'Nonzero DWT\ncoefficients\n(N={nnz_dwt})')
77
axes[2].imshow(fswavedecn_result.coeffs != 0, **imshow_kwargs)
78
axes[2].set_title(f'Nonzero FSWT\ncoefficients\n(N={nnz_fswavedecn})')
79
for ax in axes:
80
    ax.set_axis_off()
81

82
plt.show()
83