Wavelet and DCT Image Compression

Compare wavelet and discrete cosine representations for image approximation. Reconstruct images at different coefficient budgets to see how localization, block structure, and thresholding influence compression artifacts.

Run this tour

Run the cells in order with a Python 3 kernel. The first cell locates the companion data and toolbox and installs missing dependencies when needed. All worked examples include their implementation directly in this notebook. Random seeds make comparisons reproducible; you can change them to explore other samples.

# Locate the companion toolbox locally, or fetch it for a standalone/Colab copy.
from pathlib import Path
import importlib.util
import os
import subprocess
import sys

working = Path.cwd()
candidates = [working, working / "python", working.parent / "python"]
python_dir = next((p for p in candidates if (p / "nt_toolbox").is_dir()), None)
if python_dir is None:
    checkout = working / "numerical-tours-support"
    if not checkout.exists():
        subprocess.run(
            [
                "git",
                "clone",
                "--depth",
                "1",
                "--branch",
                "master",
                "https://github.com/gpeyre/numerical-tours.git",
                str(checkout),
            ],
            check=True,
        )
    python_dir = checkout / "python"
os.chdir(python_dir)
if str(python_dir) not in sys.path:
    sys.path.insert(0, str(python_dir))
requirements = python_dir / "requirements.txt"
if any(
    importlib.util.find_spec(name) is None
    for name in [
        "numpy",
        "scipy",
        "matplotlib",
        "skimage",
        "sklearn",
        "pywt",
        "ipywidgets",
        "cvxpy",
        "skfmm",
        "autograd",
        "progressbar",
        "celer",
    ]
):
    subprocess.run(
        [sys.executable, "-m", "pip", "install", "-r", str(requirements)], check=True
    )

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(0)
plt.rcParams.update(
    {
        "figure.figsize": (8, 4),
        "figure.dpi": 100,
        "axes.spines.top": False,
        "axes.spines.right": False,
        "font.size": 11,
        "image.cmap": "gray",
    }
)
%matplotlib inline

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");

References and further reading

• Stéphane Mallat. A Wavelet Tour of Signal Processing: The Sparse Way. 2009, 3rd ed., Academic Press. Multiresolution analysis, sparse approximation, and wavelet algorithms.

• Ingrid Daubechies. Ten Lectures on Wavelets. 1992, SIAM. Compactly supported orthogonal wavelets and their regularity.

• Gabriel Peyré. Advanced Signal, Image and Surface Processing. 2010, course notes. A mathematical companion to the Numerical Tours.

• David L. Donoho. De-noising by Soft-Thresholding. 1995, IEEE Transactions on Information Theory 41(3), 613–627. Why shrinkage of wavelet coefficients suppresses noise.

• Pauli Virtanen et al.. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. 2020, Nature Methods 17, 261–272. The numerical routines used for transforms, interpolation, and optimization.