Numerical Tours
← All Python tours
Download notebook Open in Colab

Explore Fourier and Wavelets Interactively

Use interactive controls to vary how many transform coefficients an image keeps. Compare reconstruction quality across Fourier and wavelet representations, and connect abrupt visual changes with the coefficients being discarded.

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 notebook compares 3 basic ways to approximate (and thus in some loose sense "compress") images:

  • Reducing the number of pixels by averaging on squares,
  • Keeping the low-frequency Fourier coefficients,
  • Keeping the largest Wavelet coefficients.
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
import warnings
import pywt  # to compute wavelet transforms

Load an image from the internet. The image has range of value [0,1].

def preprocess(f, n):
    n0 = np.min(f.size)
    f = f.crop((0, 0, n0, n0))
    f = f.resize((n, n))
    f = np.array(f)
    f = np.clip(f / np.max(f.flatten()), 0, 1)
    return f

n = 512
f = preprocess(Image.open("nt_toolbox/data/flowers.png").convert("RGB"), n)

plt.rcParams["figure.figsize"] = [8, 5]
plt.imshow(f)
plt.axis("off");
No description has been provided for this image 
def approx_pixels(f, r):
    p = int(n * np.sqrt(r))
    f = Image.fromarray(np.uint8(f * 255))
    f = f.resize((p, p), Image.NEAREST)
    f = f.resize((n, n), Image.NEAREST)
    f = np.array(f) / 255.0
    return np.clip(f, 0, 1)

def approx_fourier(f, r):
    m = int(r * n * n)
    # compute a 1D indexing from low to high frequencies
    x = np.concatenate((np.arange(0, int(n / 2)), np.arange(-int(n / 2), 0)))
    D = x[:, None] ** 2 + x[None, :] ** 2
    mask = D < np.sort(D.flatten())[m]
    # inverse the masked FFT
    F = np.fft.fft2(f, axes=(0, 1))
    F = F * mask[:, :, None]
    fM = np.real(np.fft.ifft2(F, axes=(0, 1)))
    return np.clip(fM, 0, 1), np.fft.fftshift(mask)

def approx_wavelets(f, r):
    m = int(r * n * n)
    wavetype = "bior3.5"
    nlvel = int(np.floor(np.log2(n)) - 3)
    f1 = np.zeros(f.shape)
    for k in range(3):
        with warnings.catch_warnings():
            warnings.simplefilter(
                "ignore"
            )  # pywlet issue a warning for too low level decomposition
            fW = pywt.wavedec2(f[:, :, k], wavetype, level=nlvel)
        # turn coefficient into a convenient array
        fWI, coeff_slices = pywt.coeffs_to_array(fW)
        D = np.sort(np.abs(fWI.flatten()))  # sorted coefficients magnitude
        fWI = fWI * (abs(fWI) > D[-m - 1])  # threshold
        fW1 = pywt.array_to_coeffs(fWI, coeff_slices, output_format="wavedec2")
        f1[:, :, k] = pywt.waverec2(fW1, wavetype)
    mask = abs(fWI) > D[-m - 1]
    return np.clip(np.array(f1), 0, 1), mask

def render_display(r):
    fP = approx_pixels(f, r)
    fF, mF = approx_fourier(f, r)
    fW, mW = approx_wavelets(f, r)
    plt.rcParams["figure.figsize"] = [16, 10]
    plt.subplot(2, 3, 1)
    plt.imshow(fP)
    plt.axis("off")
    plt.subplot(2, 3, 2)
    plt.imshow(fF)
    plt.axis("off")
    plt.subplot(2, 3, 3)
    plt.imshow(fW)
    plt.axis("off")
    plt.subplot(2, 3, 5)
    plt.imshow(mF)
    plt.axis("off")
    plt.subplot(2, 3, 6)
    plt.imshow(mW)
    plt.axis("off")


r = 0.01  # compressio ratio
render_display(r)
No description has been provided for this image 
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets

@interact(r=(0.01 / 2, 0.1, 0.01 / 2))
def h(r=0.05):
    render_display(r)
No description has been provided for this image

References and further reading