Image Approximation in Orthogonal Bases¶
Compare linear and nonlinear approximations in several orthogonal bases. By keeping a fixed coefficient budget, we can measure how efficiently each basis represents smooth regions, textures, and edges.
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
$\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$
This numerical tour uses several orthogonal bases to perform non-linear image approximation.
import numpy as np
import scipy as scp
import pylab as pyl
import matplotlib.pyplot as plt
from nt_toolbox.general import clamp, np, plt, pylab, rescale
from nt_toolbox.signal import (
imageplot,
load_image,
np,
perform_wavortho_transf,
plot_wavelet,
plt,
pylab,
snr,
)
import warnings
%matplotlib inline
Best $M$-terms Non-linear Approximation¶
This tours makes use of an orthogonal base $ \Bb = \{ \psi_m \}_{m=0}^{N-1} $ of the space $\RR^N$ of the images with $N$ pixels.
The best $M$-term approximation of $f$ is obtained by a non-linear thresholding
$$ f_M = \sum_{ \abs{\dotp{f}{\psi_m}}>T } \dotp{f}{\psi_m} \psi_m, $$ where the value of $T>0$ should be carefully selected so that only $M$ coefficients are not thresholded, i.e.
$$ \abs{ \enscond{m}{ \abs{\dotp{f}{\psi_m}}>T } } = M. $$
The goal is to use an ortho-basis $ \Bb $ so that the error $ \norm{f-f_M} $ decays as fast as possible when $M$ increases, for a large class of images.
This tour studies several different orthogonal bases: Fourier, wavelets (which is at the heart of JPEG-2000), cosine, local cosine (which is at the heart of JPEG).
First we load an image of $ N = n \times n $ pixels.
n = 512
f = rescale(load_image("nt_toolbox/data/hibiscus.bmp", n))
Display it.
plt.figure(figsize=(5, 5))
imageplot(f)
Fourier Approximation¶
The discrete 2-D Fourier atoms are defined as: $$ \psi_m(x) = \frac{1}{\sqrt{N}} e^{ \frac{2i\pi}{n} ( x_1 m_1 + x_2 m_2 ) }, $$ where $ 0 \leq m_1,m_2 < n $ indexes the frequency.
The set of inner products $ \{ \dotp{f}{\psi_m} \}_m $ is computed in $O(N \log(N))$ operations with the 2-D Fast Fourier Transform (FFT) algorithm (the pylab function is fft2).
Compute the Fourier transform using the FFT algorithm. Note the normalization by $1/\sqrt{N}$ to ensure orthogonality (energy conservation) of the transform.
fF = pyl.fft2(f) / n
Display its magnitude (in log scale). We use the pylab function fftshift to put the low frequency in the center.
plt.figure(figsize=(5, 5))
imageplot(np.log(1e-5 + abs(pyl.fftshift(fF))))
An image is recovered from a set of coefficients $c_m$ using the inverse Fourier Transform (pylab function ifft2)) that implements the formula
$$ f_M = \sum_m c_m \psi_m. $$
Perform a thresholding.
T = 0.3
c = np.multiply(fF, (abs(fF) > T))
Inverse the Fourier transform.
fM = np.real(pyl.ifft2(c) * n)
Display the approximation.
plt.figure(figsize=(5, 5))
imageplot(clamp(fM))
Worked example 1
Compute a best $M$-term approximation in the Fourier basis of $f$, for $M \in \{N/100, N/20\}$. Compute the approximation using a well chosen hard threshold value $T$.
plt.figure(figsize=(8, 8))
Mlist = [round(n**2 / 100), round(n**2 / 20)]
for i in range(len(Mlist)):
M = Mlist[i]
# perform hard thresholding
F = pyl.fft2(f)
a = np.sort(np.ravel(abs(F)))[::-1] # sort a 1D copy of F in descending order
T = a[M]
fFT = np.multiply(F, (abs(F) > T))
fM = np.real(pyl.ifft2(fFT))
imageplot(
clamp(fM), "M/N = %.2f, SNR = %.1f dB" % (M / n**2, snr(f, fM)), [1, 2, i + 1]
)
The best $M$-term approximation error is computed using the conservation of energy as
$$ \epsilon[M]^2 = \norm{f-f_M}^2 = \sum_{ \abs{\dotp{f}{\psi_m}} \leq T } \abs{\dotp{f}{\psi_m}}^2. $$
If one denotes by $ \{ c_R[k] \}_{k=0}^{N-1} $ the set of coefficients magnitudes $ \abs{\dotp{f}{\psi_m}} $ ordered by decaying magnitudes, then this error is easily computed as $$ \epsilon[M]^2 = \sum_{k=M}^{N-1} c_R[k]^2 = \norm{f}^2 - \sum_{k=0}^{M-1} c_R[k]^2. $$ This means that $\epsilon^2$ is equal to $\norm{f}^2$ minus the discrete primitive of $ c_R^2 $.
Worked example 2
Compute and display in log scales the ordered coefficients $c_R$. Hint: a discrete primitive can be computed using the numpy function cumsum.
cR = np.sort(np.ravel(abs(fF)))[::-1]
h = plt.plot(np.log10(cR), linewidth=2)
plt.xlim(0, n**2)
plt.show()
Worked example 3
Compute and display in log-scale the non-linear approximation error $\epsilon[M]^2$. Store the values of $\epsilon[M]^2$ in a vector $err\_fft$.
from pylab import linalg
err_fft = [max(e, 1e-10) for e in linalg.norm(f) ** 2 - np.cumsum(cR**2)]
plt.plot(np.log10(err_fft / linalg.norm(f) ** 2), linewidth=2)
plt.title(r"$\log_{10}(\epsilon^2[M]/ ||f||^2)$")
plt.xlim(1, n**2 / 50)
plt.ylim(-2.35, 0)
plt.show()
Wavelet Approximation¶
The Wavelet basis of continuous 2-D functions is defined by by scaling and translating three mother atoms $ \{\psi^H,\psi^V,\psi^D\} $: $$ \psi_{j,n}^k(x) = \frac{1}{2^j}\psi^k\pa{\frac{x-2^j n}{2^j}} $$
Non-linear wavelet approximation is a the heart of the JPEG-2000 compression standard.
The set of inner products $ \{ \dotp{f}{\psi_m} \}_m $ is computed in $O(N)$ operations with the 2-D Fast Wavelet Transform algorithm.
Perform a wavelet transform. Here we use a daubechies wavelet transform.
from numpy import linalg
from nt_toolbox.compute_wavelet_filter import compute_wavelet_filter, np
Jmin = 1
h = compute_wavelet_filter("Daubechies", 10)
fW = perform_wavortho_transf(f, Jmin, +1, h)
Display the coefficients.
plt.figure(figsize=(8, 8))
plot_wavelet(fW, Jmin)
plt.title("Wavelet coefficients")
plt.show()
Worked example 4
Compute a best $M$-term approximation in the wavelet basis of $f$, for $M \in \{N/100, N/20\}$. Compute the approximation using a well chosen hard threshold value $T$. Note that the inverse wavelet transform is obtained by replacing the +1 by a -1 in the definition of the transform.
plt.figure(figsize=(8, 8))
Mlist = [round(n**2 / 100), round(n**2 / 20)]
for i in range(len(Mlist)):
M = Mlist[i]
a = np.sort(np.ravel(abs(fW)))[::-1] # sort a 1D copy of F in descending order
T = a[M]
fWT = np.multiply(fW, (abs(fW) > T))
fM = np.real(perform_wavortho_transf(fWT, Jmin, -1, h))
imageplot(
clamp(fM), "M/N = %.2f, SNR = %.1f dB" % (M / n**2, snr(f, fM)), [1, 2, i + 1]
)
Worked example 5
Compute and display in log-scale the non-linear approximation error $\epsilon[M]^2$. Compares the Fourier and wavelets approximations. Store the values of $\epsilon[M]^2$ in a vector $err\_wav$.
from pylab import linalg
cR = np.sort(np.ravel(abs(fW)))[::-1]
err_wav = [max(e, 1e-10) for e in linalg.norm(f) ** 2 - np.cumsum(cR**2)]
plt.plot(
np.log10(err_fft / linalg.norm(f) ** 2), linewidth=2, color="red", label="Fourier"
)
plt.plot(
np.log10(err_wav / linalg.norm(f) ** 2), linewidth=2, color="blue", label="Wavelets"
)
plt.title(r"$\log_{10}(\epsilon^2[M]/ ||f||^2)$")
plt.xlim(1, n**2 / 50)
plt.ylim(-2.35, 0)
plt.legend()
plt.show()
Cosine Approximation¶
The discrete cosine approximation (DCT) is similar to the Fourier approximation, excepted that it used symmetric boundary condition instead of periodic boundary condition, and is thus more useful to approximate image.
A 1-D cosine atom of $n$ sample is defined as $$ \bar\psi_m(x) = \frac{1}{\sqrt{N}} \cos\pa{ \frac{2\pi}{N} (x-1/2) m } $$ A 2-D cosine atom is obtained by tensor product of 1-D atoms $$ \psi_{m_1,m_2}(x_1,x_2) = \bar\psi_{m_1}(x_1) \bar\psi_{m_2}(x_2). $$ On the contrary to the Fourier 2-D atoms, these 2-D DCT atoms are not oriented (they contains 4 Fourier frequencies).
The set of inner products $ \{ \dotp{f}{\psi_m} \}_m $ is computed in $O(N \log(N))$ operations with the 2-D Fast Cosine Transform algorithm.
from scipy import fftpack
def dct2(f):
return np.transpose(
fftpack.dct(np.transpose(fftpack.dct(f, norm="ortho")), norm="ortho")
)
def idct2(f):
return np.transpose(
fftpack.idct(np.transpose(fftpack.idct(f, norm="ortho")), norm="ortho")
)
fC = dct2(f)
Display the magnitude of the DCT coefficients. Note that the low frequencies are in the upper-left corner.
plt.figure(figsize=(5, 5))
imageplot(np.log(1e-5 + abs(fC)))
Worked example 6
Compute a best $M$-term approximation in the wavelet basis of $f$, for $M \in \{N/100, N/20\}$. Compute the approximation using a well chosen hard threshold value $T$. Note that the inverse DCT transform is obtained with the function idct.
plt.figure(figsize=(8, 8))
Mlist = [round(n**2 / 100), round(n**2 / 20)]
for i in range(len(Mlist)):
M = Mlist[i]
# perform hard thresholding
a = np.sort(np.ravel(abs(fC)))[::-1] # sort a 1D copy of F in descending order
T = a[M]
fCT = np.multiply(fC, abs(fC) > T)
fM = idct2(fCT)
imageplot(
clamp(fM), "M/ N = %.2f, SNR = %.1f dB" % (M / n**2, snr(f, fM)), [1, 2, i + 1]
)
Worked example 7
Compute and display in log-scale the non-linear approximation error $\epsilon[M]^2$. Compares the Fourier and DCT approximations. Store the values of $\epsilon[M]^2$ in a vector $err\_dct$.
from pylab import linalg
cR = np.sort(np.ravel(abs(fC)))[::-1]
err_dct = [max(e, 1e-10) for e in linalg.norm(f) ** 2 - np.cumsum(cR**2)]
plt.plot(
np.log10(err_fft / linalg.norm(f) ** 2), linewidth=2, color="red", label="Fourier"
)
plt.plot(
np.log10(err_wav / linalg.norm(f) ** 2), linewidth=2, color="green", label="DCT"
)
plt.title(r"$\log_{10}(\epsilon^2[M]/ ||f||^2)$")
plt.xlim(1, n**2 / 50)
plt.ylim(-2.35, 0)
plt.legend()
plt.show()
Local Cosine Approximation¶
To improve the global DCT approximation, one can approximate independantly small patches in the image. This corresponds to a decomposition in a local cosine basis, which is at the heart of the JPEG image compression standard.
The only parameter of the transform is the size of the square.
w = 16
Initialize at zero the transformed image in the local DCT basis.
fL = np.zeros([n, n])
Example of patch index.
i = 5
j = 7
For a given path index $(i,j)$, we extract a $(w,w)$ patch.
P = f[(i - 1) * w : i * w, (j - 1) * w : j * w]
Compute the Cosine transform of the patch using the fast DCT algorithm.
fL[(i - 1) * w : i * w, (j - 1) * w : j * w] = dct2(P)
Display the patch and its coefficients. We removed the low frequency of $P$ for display purpose only.
plt.figure(figsize=(8, 8))
imageplot(P, "Patch", [1, 2, 1])
imageplot(dct2(P - np.mean(P)), "DCT", [1, 2, 2])
Worked example 8
Compute the local DCT transform $f_L$ by transforming each patch.
fL = np.zeros([n, n])
for i in range(1, n // w + 1):
for j in range(1, n // w + 1):
fL[(i - 1) * w : i * w, (j - 1) * w : j * w] = dct2(
f[(i - 1) * w : i * w, (j - 1) * w : j * w]
)
Display the coefficients.
plt.figure(figsize=(5, 5))
imageplot(np.clip(abs(fL), 0, 0.005 * w * w))
Worked example 9
Compute the inverse local DCT transform of the coefficients $f_L$ by inverse transforming each patch using the function idct2.
from numpy import linalg
f1 = np.copy(fL)
for i in range(1, n // w + 1):
for j in range(1, n // w + 1):
f1[(i - 1) * w : i * w, (j - 1) * w : j * w] = idct2(
f1[(i - 1) * w : i * w, (j - 1) * w : j * w]
)
print("Error |f-f1|/|f| =", linalg.norm(f - f1) / linalg.norm(f))
Worked example 10
Compute a few best $M$-term approximations in the Local DCT basis of $f$.
plt.figure(figsize=(8, 8))
for u in range(len(Mlist)):
M = Mlist[u]
a = np.sort(np.ravel(abs(fL)))[::-1] # sort a 1D copy of F in descending order
T = a[M]
fLT = np.multiply(fL, abs(fL) > T)
# fLT = perform_thresholding(fL, M, 'largest')
fM = fLT
for i in range(1, n // w + 1):
for j in range(1, n // w + 1):
fM[(i - 1) * w : i * w, (j - 1) * w : j * w] = idct2(
fLT[(i - 1) * w : i * w, (j - 1) * w : j * w]
)
imageplot(
clamp(fM), "M/N = %.2f, SNR = %.1f dB" % (M / n**2, snr(f, fM)), [1, 2, u + 1]
)
Worked example 11
Compute and display in log-scale the non-linear approximation error $\epsilon[M]^2$. Store the values of $\epsilon[M]^2$ in a vector |err_ldct|. Compares the Fourier, Wavelets, DCT and local-DCT approximations.
from pylab import linalg
cR = np.sort(np.ravel(abs(fL)))[::-1]
err_ldct = [max(e, 1e-10) for e in linalg.norm(f) ** 2 - np.cumsum(cR**2)]
plt.plot(
np.log10(err_fft / linalg.norm(f) ** 2), linewidth=2, color="red", label="Fourier"
)
plt.plot(
np.log10(err_wav / linalg.norm(f) ** 2), linewidth=2, color="blue", label="Wavelets"
)
plt.plot(
np.log10(err_dct / linalg.norm(f) ** 2), linewidth=2, color="purple", label="DCT"
)
plt.plot(
np.log10(err_ldct / linalg.norm(f) ** 2),
linewidth=2,
color="orange",
label="Local DCT",
)
plt.title(r"$\log_{10}(\epsilon^2[M]/ ||f||^2)$")
plt.xlim(1, n**2 / 50)
plt.ylim(-2.35, 0)
plt.legend()
plt.show()
Comparison of Wavelet Approximations of Several Images¶
An image is more complicated than an other one for a given orthogonal basis if its approximation error decays more slowly.
First load several high resolution images.
n = 512
fList = np.zeros([n, n, 4])
fList[:, :, 0] = rescale(load_image("nt_toolbox/data/regular3.bmp", n))
fList[:, :, 1] = rescale(load_image("nt_toolbox/data/phantom.bmp", n))
fList[:, :, 2] = rescale(load_image("nt_toolbox/data/lena.bmp", n))
fList[:, :, 3] = rescale(load_image("nt_toolbox/data/mandrill.bmp", n))
Display them.
plt.figure(figsize=(7, 7))
for i in range(4):
imageplot(fList[:, :, i], "", [2, 2, i + 1])
Worked example 12
Compare the approximation error decay for those images. Display $ \log_{10}(\norm{f-f_M}) $ as a function of $\log_{10}(M)$.
plt.figure(figsize=(8, 5))
names = ["regular3", "phantom", "lena", "mandrill"]
for i in range(np.shape(fList)[2]):
fW = perform_wavortho_transf(fList[:, :, i], Jmin, +1, h)
cR = np.sort(np.ravel(abs(fW)))[::-1]
err = [e for e in linalg.norm(fList[:, :, i]) ** 2 - np.cumsum(cR**2)]
Err = err[10 : n * n // 10]
plt.plot(
np.log10(np.arange(10, n * n // 10)),
np.log10(Err / Err[0]),
linewidth=2,
label=names[i],
)
plt.title(r"$\log_{10}(\epsilon^2[M]/ ||f||^2)$")
plt.xlim(1, np.log10(n * n // 10))
plt.ylim(-7, 0)
plt.legend(loc=3)
plt.show()
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.