This numerical tour showcases a comparison of Wavelet and DCT image "compression" by thresholding the transformed coefficients. This corresponds to the best non-linear approximation in these orthogonal bases.
import numpy as np
import scipy
import matplotlib.pyplot as plt
import pywt
import pywt.data
from scipy.fftpack import dct, idct
Load and display the original image.
f = pywt.data.camera()
n = f.shape[0] # size of the image
plt.imshow(f, cmap='gray')
plt.axis('off');
Define the thresholding operator. Keep the $M$ largest coefficients of a vector.
def thresh(fw,M):
a = np.sort(np.ravel(abs(fw)))[::-1] #sort a 1D copy of fw in descending order
T = a[M]
return fw * (abs(fw)>T)
Define the 2D DCT transform by applying the 1D DCT along each axis.
def dct2(a):
return dct( dct( a, axis=0, norm='ortho' ), axis=1, norm='ortho' )
def idct2(a):
return idct( idct( a, axis=1 , norm='ortho'), axis=0 , norm='ortho')
Usual measure of quality (of course debatable) is the PSNR (log of the $\ell^2$ error).
def psnr(img1, img2):
mse = np.mean( (img1 - img2) ** 2 )
PIXEL_MAX = 255.0
return 20 * np.log10(PIXEL_MAX / np.sqrt(mse))
This controls the compression ration. Set ratio=100/100 to get the full image.
ratio = 1.5/100
M = int(ratio * n**2)
Perform the DCT compression.
fd = dct2(f);
fd1 = thresh(fd,M)
f1 = idct2(fd1);
Type of wavelet transform -- you can change this to see the effect
# WT = 'db2'; # celebrated Daubechies Wavelet
# WT = 'bior3.5'; # biorthogonal but maybe a bit too wide
# WT = 'haar'; # produces a blocky approximation
WT = 'bior2.4'; # not exactly orthogonal but considered good
Performs the Wavelet compression.
fw, S = pywt.coeffs_to_array( pywt.wavedec2(f, WT) );
fw1 = thresh(fw,M)
f2 = pywt.waverec2( pywt.array_to_coeffs(fw1,S,output_format='wavedec2'), WT)
Display DCT and Wavelet compression.
plt.subplot(1, 2, 1)
plt.imshow(fd1!=0); plt.axis('off'); plt.title('DCT coefs');
plt.subplot(1, 2, 2)
plt.imshow(np.clip(f1,0,255), cmap='gray'); plt.axis('off');
plt.title('Compres. PSNR=' + '{:.2f}'.format(psnr(f,f1)) + 'dB' );
plt.subplot(1, 2, 1)
plt.imshow(fw1!=0); plt.axis('off'); plt.title('Wavelets coefs');
plt.subplot(1, 2, 2)
plt.imshow(np.clip(f2,0,255), cmap='gray'); plt.axis('off');
plt.title('Compres. PSNR=' + '{:.2f}'.format(psnr(f,f2)) + 'dB' );
plt.imshow(np.clip(f1,0,255), cmap='gray'); plt.axis('off'); plt.title('DCT.compression');
plt.imshow(np.clip(f2,0,255), cmap='gray'); plt.axis('off'); plt.title('Wav.compression');