In the previous chapters, all operations were performed in the spatial domain, where algorithms act directly on pixel intensity values.
This chapter presents a complementary approach: the frequency domain, in which the image is represented by the spatial variations of intensity, rather than only by the individual pixel values.
The concept of spatial frequency describes how quickly intensity varies across the image. Slow variations correspond to low frequencies, while edges, fine details, and noise correspond to high frequencies.
This representation relies on the fact that any discrete digital image can be decomposed into a combination of orthogonal functions. The Fourier Transform uses a basis of two-dimensional complex exponentials (equivalent to sinusoids with specific orientation and frequency). Other transforms, such as the Cosine Transform (DCT) and the Wavelet Transform (DWT), use different families of basis functions — two-dimensional cosines in the case of the DCT, and functions with compact support in the case of wavelets.
Among the main applications of this representation are:
By the end of this chapter, you will be able to:
import os, urllib.request
url = "https://raw.githubusercontent.com/fzampirolli/pdi-vc/master/morph/config.py"
if not os.path.exists("config.py"):
urllib.request.urlretrieve(url, "config.py")
import config
config.setup()
from morph import mm
import numpy as np✅ Environment ready. Morph: 1.1.9 | OpenCV: 5.0.0Fourier analysis is based on the principle that any periodic signal can be represented as a sum of sinusoidal functions with different frequencies, amplitudes, and phases. This concept also applies to digital images, allowing them to be represented in the frequency domain instead of the spatial domain.
Figure 5.1 illustrates this decomposition for a one-dimensional signal. In the case of an image, the Discrete Fourier Transform (DFT) converts the intensity matrix \(f(x,y)\) into a set of coefficients that describes the contribution of the different spatial frequencies present in the image.
Before studying two-dimensional images, the simulator in Figure 5.2 illustrates the principle of Fourier analysis for one-dimensional signals: a waveform can be approximated by the sum of sinusoids with different frequencies and amplitudes.
As new terms are added, the sum of the sinusoids (black curve) approaches the reference waveform (dashed). The lower plot displays the amplitude spectrum, indicating the contribution of each frequency to the reconstruction of the signal.
Explore the simulator and answer:
1. How many terms are needed for a good approximation of the square wave?
With approximately 15 to 20 terms, the waveform already closely approximates the reference. However, near the discontinuities, a small oscillation remains, known as the Gibbs phenomenon, which does not disappear even with the addition of more terms.
2. Which waveform converges most rapidly? Why?
The triangular wave converges most rapidly, because the amplitudes of its harmonics decay faster than those of the square wave and the sawtooth wave. As a result, few terms already produce a good approximation.
3. How does the spectrum change between the square wave and the triangular wave?
Both contain only odd harmonics, but in the triangular wave, the amplitudes decrease much more rapidly. Thus, few harmonics are sufficient to reconstruct the signal with good accuracy.
When applying the Discrete Fourier Transform (DFT) to an image and visualizing the magnitude of its coefficients (see Figure 5.5), one obtains the magnitude spectrum, which shows the distribution of spatial frequencies present in the image.
The coefficient located at the origin of the DFT, known as the DC component (Direct Current), corresponds to the zero frequency and represents the average intensity of the image. By convention, this coefficient is stored in the upper left corner of the spectrum. To facilitate its interpretation, the FFT Shift operation is applied, which shifts the DC component to the center of the image. After this shift, low frequencies are concentrated in the central region, while high frequencies are near the edges, as summarized in Table 5.1.
| Spectrum region | Predominant components | Examples in the image |
|---|---|---|
| Center (low frequencies) | Slow spatial variations | Illumination, homogeneous regions, and global shapes |
| Intermediate region (mid frequencies) | Intermediate-scale variations | Textures and repetitive patterns |
| Edges (high frequencies) | Rapid spatial variations | Contours, fine details, and noise |
This organization facilitates the interpretation of the spectrum and the design of filters. Attenuating low frequencies reduces global intensity variations, while attenuating high frequencies smooths the image by reducing fine details and part of the noise.
Before presenting the mathematical formulation of the Discrete Fourier Transform (DFT), it is useful to analyze its inverse, called the Inverse Discrete Fourier Transform (IDFT). Consider a spectrum in which all coefficients are zero, except one. An example of this construction is presented in the code of Figure 5.3 and can be explored interactively in the simulator of Figure 5.4.
The reconstructed image is a two-dimensional sinusoid. The position of the coefficient in the spectrum determines its orientation and its spatial frequency, while its magnitude and phase define, respectively, its amplitude and spatial displacement. Thus, each DFT coefficient represents a sinusoidal component, and the original image can be reconstructed by summing all these components.
import cv2
N_grid = 100
espectro_vazio = np.zeros((N_grid, N_grid), dtype=complex)
# Turning on a single point (frequency) off-center
u0, v0 = 10, 5
espectro_vazio[N_grid//2 - v0, N_grid//2 - u0] = 1000
# Returning to the spatial domain (IDFT)
onda_2d = np.real(np.fft.ifft2(np.fft.ifftshift(espectro_vazio)))
onda_vis = cv2.normalize(onda_2d, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)
espectro_vis = cv2.normalize(
np.abs(espectro_vazio), None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)
# Visual highlight of the point
espectro_color = cv2.cvtColor(espectro_vis, cv2.COLOR_GRAY2BGR)
cv2.circle(espectro_color, (N_grid//2 - u0, N_grid//2 - v0), 2, (0, 0, 255), -1)
mm.show([espectro_color, onda_vis],
titles=["Spectrum (1 active point)", "Resulting 2D Wave (IDFT)"],
cols=2, figsize=(10, 4))
Consider an image \(f(x,y)\) with dimensions \(M \times N\). Its 2D Discrete Fourier Transform (DFT) is defined by:
\[ F(u,v) = \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x,y)\, e^{-j2\pi\left(\frac{ux}{M}+\frac{vy}{N}\right)} \tag{5.1}\]
where \(u = 0, 1, \ldots, M-1\) and \(v = 0, 1, \ldots, N-1\) represent the discrete frequencies in the horizontal and vertical directions, respectively. The exponential term corresponds to a two-dimensional sinusoid, whose frequency and orientation are determined by the indices \((u,v)\).
The 2D Inverse Discrete Fourier Transform (IDFT) reconstructs the original image from its coefficients:
\[ f(x,y) = \frac{1}{MN} \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} F(u,v)\, e^{j2\pi\left(\frac{ux}{M}+\frac{vy}{N}\right)} \tag{5.2}\]
Equations Equation 5.1 and Equation 5.2 show that the DFT and the IDFT form a pair of transformations: the former converts the image into the frequency domain, while the latter exactly reconstructs the original image from its coefficients.
The term \(j\) denotes the imaginary unit, defined by \(j^2 = -1\). In engineering and signal processing, \(j\) is adopted instead of \(i\) to avoid conflict with the notation for electric current. Its use in the complex exponential, governed by Euler’s formula (\(e^{j\theta} = \cos\theta + j\sin\theta\)), allows for a compact representation of the amplitude and phase of each spatial frequency present in the image.
The coefficient \(F(0,0)\), termed the DC component (Direct Current), equals the sum of the intensities of all pixels in the image (see Figure 5.5):
\[ F(0,0)=MN\,\bar{f}, \]
where \(\bar{f}\) is the average intensity of the image. Therefore, the DC component represents the mean intensity level and, in most natural images, holds the largest magnitude of the spectrum.
The remaining coefficients represent variations around this mean. After applying the FFT Shift, the DC component is moved to the center of the spectrum, concentrating low frequencies in the central region and high frequencies at the edges.
Each coefficient of the Discrete Fourier Transform (DFT) is a complex number and can be written as
\[ F(u,v)=R(u,v)+j\,I(u,v), \]
where \(R(u,v)\) and \(I(u,v)\) correspond, respectively, to the real and imaginary parts of the coefficient. From Equation Equation 5.1, we obtain
\[ R(u,v)= \sum_{x=0}^{M-1}\sum_{y=0}^{N-1} f(x,y) \cos\!\left( 2\pi\left(\frac{ux}{M}+\frac{vy}{N}\right) \right), \]
and
\[ I(u,v)= - \sum_{x=0}^{M-1}\sum_{y=0}^{N-1} f(x,y) \sin\!\left( 2\pi\left(\frac{ux}{M}+\frac{vy}{N}\right) \right). \]
From this representation, two fundamental quantities are defined:
\[ |F(u,v)|=\sqrt{R(u,v)^2+I(u,v)^2}; \]
\[ \phi(u,v)=\operatorname{atan2}\!\left(I(u,v),\,R(u,v)\right). \]
Thus, each coefficient can also be written in its polar form,
\[ F(u,v)=|F(u,v)|\,e^{j\phi(u,v)}. \]
The Fourier spectrum can, therefore, be visualized by means of two distinct images: the magnitude spectrum, typically used to analyze the distribution of frequencies, and the phase spectrum, which describes the spatial organization of the sinusoidal components.
Although the magnitude spectrum is the most used for visual inspection, the phase contains a large portion of the structural information of the image. The combination of magnitude and phase allows the exact reconstruction of the original image via the IDFT.
A classic demonstration consists of combining the magnitude of one image with the phase of another and reconstructing the result. This experiment highlights that:
When an image is reconstructed with the magnitude of A and the phase of B, the result tends to resemble B more than A, showing that phase is the main component responsible for the spatial organization of the scene. However, magnitude remains important, as it modulates the contrast of the reconstructed structures. Thus, a faithful reconstruction depends on the consistent combination of magnitude and phase.
An example of this behavior is presented in Figure 5.6.
The phase of a signal plays distinct roles in audio and images:
In both domains, magnitude is related to the intensity of frequency components: in audio, it influences timbre and perceived loudness; in images, it influences contrast and texture.
# ── Experiment: The Importance of Phase ───────────────────────────────────────
# ── Image Loading ────────────────────────────────────────────────────────────
url = "https://upload.wikimedia.org/wikipedia/commons/2/25/GAZI.MD.AHAD_11.jpg"
caminho = "imagens/coins.jpg"
if not os.path.exists(caminho):
os.makedirs("imagens", exist_ok=True)
img_obj = mm.read(url, pil=True)
mm.write(img_obj, caminho)
else:
img_obj = mm.read(caminho, pil=True)
img_color = np.array(img_obj)
img_gray = mm.gray(img_color)
img_a = cv2.resize(img_gray, (400, 400))
# Create a synthetic B image (geometric pattern)
img_b = np.zeros((400, 400), dtype=np.uint8)
cv2.rectangle(img_b, (100, 100), (300, 300), 255, -1)
cv2.circle(img_b, (200, 200), 150, 128, 10)
FA = np.fft.fft2(img_a)
FB = np.fft.fft2(img_b)
# Phase Swap
rec_A_mag_B_fase = np.real(np.fft.ifft2(np.abs(FA) * np.exp(1j * np.angle(FB))))
rec_B_mag_A_fase = np.real(np.fft.ifft2(np.abs(FB) * np.exp(1j * np.angle(FA))))
mm.show(
[img_a, img_b, rec_A_mag_B_fase, rec_B_mag_A_fase],
titles=["Image A", "Image B", "Mag(A) + Phase(B)", "Mag(B) + Phase(A)"],
cols=4, figsize=(16, 4)
)
print("💡 Phase preserves edges and contours; magnitude controls contrast and")
print("texture. In stereo audio, phase affects spatial location; in")
print("images, it determines the organization of the scene.")
💡 Phase preserves edges and contours; magnitude controls contrast and
texture. In stereo audio, phase affects spatial location; in
images, it determines the organization of the scene.The Convolution Theorem establishes a fundamental relationship between the spatial and frequency domains:
\[ f(x,y) \circledast h(x,y) \;\overset{\mathcal{F}}{\longleftrightarrow}\; F(u,v)\,H(u,v) \tag{5.3}\]
where \(\circledast\) denotes discrete circular convolution. Thus, the convolution between an image \(f(x,y)\) and a filter \(h(x,y)\) can be replaced by the multiplication of their spectra.
In practice, to obtain the same result as the linear convolution performed in the spatial domain, zero-padding is applied before the Fast Fourier Transform (FFT), avoiding artifacts at the image borders.
However, frequency-domain filtering is not always the most efficient alternative. For filters such as the Gaussian and the Box Filter, the separability property allows a significant reduction in the computational cost of convolution in the spatial domain.
A separable kernel can be written as the outer product of two one-dimensional vectors,
\[ H = v\,h^T, \]
allowing the two-dimensional convolution to be replaced by two consecutive one-dimensional convolutions: one in the horizontal direction and another in the vertical direction.
In contrast, a non-separable kernel does not admit such a decomposition, and therefore its convolution must be performed directly over the two-dimensional neighborhood.
In practice, for a kernel of dimension \(K \times K\), direct convolution requires \(K^2\) multiplications per pixel, whereas a separable kernel requires only \(2K\) multiplications, significantly reducing the computational cost.
Consider an image of dimensions \(M \times N\) and a square filter of size \(K \times K\). Table 5.2 compares the complexity of the main filtering strategies.
| Filtering Method | Asymptotic Complexity | Dependence on \(K\) | Typical Application |
|---|---|---|---|
| Non-separable spatial | \(\mathcal{O}(MNK^2)\) | Quadratic | Small, non-separable kernels |
| Separable spatial | \(\mathcal{O}(MNK)\) | Linear | Gaussian and mean filters |
| Via FFT | \(\mathcal{O}(MN\log(MN))\) | Independent of \(K\) | Large kernels |
For small kernels, spatial convolution, especially when the filter is separable, tends to be more efficient due to the low cost of operations. As the kernel size increases, FFT-based filtering becomes more advantageous, since its cost is practically independent of the filter dimension.
The graph obtained in the test with the coin image (\(2560 \times 1920\)), presented in Figure 5.7, confirms the behavior predicted by the computational complexity analysis.
Non-separable convolution (\(\mathcal{O}(MNK^2)\))
Direct convolution exhibits quadratic growth with kernel size. For small values of \(K\), the cost is low, but it increases rapidly as the kernel grows, becoming unfeasible for real-time applications.
FFT-based filtering (\(\mathcal{O}(MN \log(MN))\))
The cost of the FFT depends only on the image size, being independent of \(K\). Therefore, its performance remains approximately constant as the kernel varies, making it advantageous for large or non-separable filters.
Separable convolution (\(\mathcal{O}(MNK)\))
Decomposing the kernel into two one-dimensional filters significantly reduces the computational cost. In practice, this approach tends to be the most efficient for separable filters, especially in optimized implementations.
In general, the choice of method depends on the size and structure of the kernel. Separable filters are more efficient in the spatial domain, while the FFT becomes more advantageous for large kernels or multiple convolutions in the frequency domain.
\[ g = \mathcal{F}^{-1}\bigl[\mathcal{F}(f)\cdot \mathcal{F}(h)\bigr] \quad \text{(FFT)} \qquad g = f \circledast h \quad \text{(direct convolution)} \qquad g = (f \circledast v) \circledast h^T \quad \text{(separable)} \tag{5.4}\]
where:
The Discrete Fourier Transform (DFT) assumes that the image is periodically extended in space, meaning that its borders repeat indefinitely.
Under this condition, multiplication in the frequency domain corresponds to a circular convolution in the spatial domain. Consequently, opposite regions of the image (top and bottom, left and right) begin to interact artificially, as illustrated in Figure 5.8.
Applying zero-padding before the FFT reduces this effect by extending the image with null values at the borders, approximating the result of linear convolution. This behavior can be interpreted in light of the Convolution Theorem, presented in Figure 5.9.
# ── Define M and N ─────────────────────────────────────────────────────────────
M, N = img_gray.shape # added line
# Simulation of a brutal shift filter
H_shift = np.zeros_like(img_gray, dtype=complex)
for u in range(M):
for v in range(N):
H_shift[u, v] = np.exp(-1j * 2 * np.pi * (u*120/M + v*120/N))
# Filtering WITHOUT padding (causes wrap-around)
F_img = np.fft.fft2(img_gray)
img_vazada = np.real(np.fft.ifft2(F_img * H_shift))
img_vazada_vis = cv2.normalize(img_vazada, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)
mm.show([img_gray, img_vazada_vis],
titles=["Original", "Filtering w/o Padding (Leakage)"], cols=2, figsize=(10, 4))
# ── Gaussian Kernel 11×11
sigma = 3.0
K = 11
ks = np.arange(K) - K // 2
gauss1d = np.exp(-ks**2 / (2 * sigma**2))
gauss1d /= gauss1d.sum()
kernel = np.outer(gauss1d, gauss1d) # separable 2D kernel
# ── Method 1: Direct spatial convolution ─────────────────────────────────────
f_float = img_gray.astype(np.float64)
conv_esp = cv2.filter2D(f_float, -1, kernel, borderType=cv2.BORDER_CONSTANT)
# ── Method 2: Frequency multiplication (via FFT) ──────────────────────────
M, N = f_float.shape
# Padding for linear convolution (avoids circular aliasing)
Mpad = 2 ** int(np.ceil(np.log2(M + K - 1)))
Npad = 2 ** int(np.ceil(np.log2(N + K - 1)))
# Positions the kernel with the origin at (0,0) and zero padding
kernel_pad = np.zeros((Mpad, Npad))
kh, kw = kernel.shape
kernel_pad[:kh, :kw] = kernel
F_img = np.fft.fft2(f_float, (Mpad, Npad))
F_kern = np.fft.fft2(kernel_pad)
conv_freq = np.real(np.fft.ifft2(F_img * F_kern))
# Crop to compensate for the shift introduced by kernel positioning
offset = K // 2
conv_freq_crop = conv_freq[offset:offset+M, offset:offset+N]
# ── Numerical verification ──────────────────────────────────────────────────
diff = np.abs(conv_esp - conv_freq_crop)
print(f"Maximum difference (|conv_spatial - conv_freq|): {diff.max():.2e}")
print(f"Mean difference (|conv_spatial - conv_freq|): {diff.mean():.2e}")
print(f"→ Convolution Theorem numerically verified.")
conv_esp_vis = cv2.normalize(conv_esp, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)
conv_freq_vis = cv2.normalize(conv_freq_crop, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)
diff_vis = cv2.normalize(diff, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)
mm.show(
[img_gray, conv_esp_vis, conv_freq_vis, diff_vis],
titles=[
"Original",
"Spatial convolution",
"Frequency multiplication",
f"Difference (max={diff.max():.1e})"
],
cols=4, figsize=(16, 5)
)Maximum difference (|conv_spatial - conv_freq|): 2.56e-13
Mean difference (|conv_spatial - conv_freq|): 2.76e-14
→ Convolution Theorem numerically verified.
The residual difference on the order of \(10^{-13}\) does not violate the Convolution Theorem, but rather reflects computational limitations inherent to floating-point arithmetic (double precision, ~\(10^{-16}\)) and to the order of operations between the two methods:
Therefore, the theoretical equality is exact, but the numerical implementation produces a practically negligible difference (relative error < \(10^{-12}\)), confirming the theorem within machine precision.
A frequency domain filter can be interpreted as a transfer function applied to the image spectrum. In this representation, each frequency coefficient is multiplied by a value between 0 and 1, which determines its attenuation or preservation. The shape of this function defines the visual effect of the filter.
Abrupt cutoff and ringing. Ideal filters with an instantaneous transition at a cutoff frequency \(D_0\) produce discontinuities in the frequency domain. This discontinuity is reflected in the spatial domain as oscillations near edges, known as ringing. This effect is associated with convolution with functions of infinite support in space, such as the sinc function, as illustrated in Figure 5.10.
Filters with smooth transition. Alternatives such as the Gaussian and Butterworth filters smooth the transition between pass and reject regions, reducing ringing. In contrast, this smoothing implies a less defined separation boundary between preserved and attenuated frequencies.
# Simulating the Ideal Frequency Filter (Cylinder) and its Spatial Representation (Sinc)
N_grid = 2**7
u = np.arange(-N_grid//2, N_grid//2)
U, V = np.meshgrid(u, u)
D = np.sqrt(U**2 + V**2)
# Frequency: Ideal Cylinder (1 at the center, 0 outside radius 20)
H_freq = np.zeros((N_grid, N_grid))
H_freq[D <= 20] = 1
# Space: The inverse results in the infamous 2D Sinc
h_space = np.fft.fftshift(np.real(np.fft.ifft2(np.fft.ifftshift(H_freq))))
fig, ax = plt.subplots(1, 2, subplot_kw={'projection': '3d'}, figsize=(12, 4))
ax[0].plot_surface(U, V, H_freq, cmap='viridis', edgecolor='none')
ax[0].set_title("Frequência: Filtro Ideal (Cilindro)")
ax[0].set_zlim(0, 1.2)
ax[1].plot_surface(U, V, h_space, cmap='plasma', edgecolor='none')
ax[1].set_title("Espaço Real: Ondulações da Sinc (Causa do Ringing)")
plt.tight_layout(); plt.show()
Low-pass filters attenuate high-frequency components, resulting in image smoothing and noise reduction. After spectrum centering (FFT Shift), the distance from each point to the center is given by:
\[ D(u,v) = \sqrt{\left(u - \tfrac{M}{2}\right)^2 + \left(v - \tfrac{N}{2}\right)^2} \tag{5.5}\]
Ideal Filter (LPFI): \[ H_{\text{ideal}}(u,v) = \begin{cases} 1, & D(u,v) \leq D_0 \\ 0, & D(u,v) > D_0 \end{cases} \tag{5.6}\]
The abrupt cutoff at \(D_0\) introduces discontinuities in the frequency domain, resulting in oscillations in the spatial domain known as ringing. This effect is associated with convolution with functions of infinite support.
Gaussian Filter (LPFG): \[ H_{\text{gauss}}(u,v) = e^{-D^2(u,v)/(2\sigma^2)} \tag{5.7}\]
The smoothness of the Gaussian function in the frequency domain avoids discontinuities, which eliminates ringing and produces a gradual transition between preserved and attenuated frequencies.
Butterworth Filter (LPFB) of order \(n\): \[ H_{\text{BW}}(u,v) = \frac{1}{1 + \left[D(u,v)/D_0\right]^{2n}} \tag{5.8}\]
The parameter \(n\) controls the smoothness of the transition between frequency pass and rejection. Small values produce smooth transitions, while large values approximate the behavior of the ideal filter, with a higher risk of ringing. A comparative example is presented in Figure 5.11.
High-pass filters can be obtained from a complementary low-pass filter, defined as:
\[ H_{\text{HP}}(u,v) = 1 - H_{\text{LP}}(u,v) \]
This type of filter preserves high-frequency components, enhancing edges and details, while attenuating regions of smooth variation.
Band-pass filters preserve only an intermediate range of frequencies, bounded by two radii \(D_L\) and \(D_H\):
\[ H_{\text{BP}}(u,v) = H_{\text{LP}}^{(D_H)}(u,v)\cdot \left[1 - H_{\text{LP}}^{(D_L)}(u,v)\right] \]
This type of filtering is useful when one wishes to simultaneously remove low- and high-frequency components, preserving only structures of intermediate scale.
An important application is the removal of periodic noise, in which regular patterns appear as localized peaks in the magnitude spectrum. These peaks can be attenuated by means of notch (band-reject) filters, positioned specifically at the undesired frequencies.
Examples of filters in the frequency domain are presented in the simulator of Figure 5.12, Figure 5.13, and Figure 5.14.
import io
def distancia_centro(M, N):
"""Distance matrix from the center of the spectrum."""
u = np.arange(M) - M // 2
v = np.arange(N) - N // 2
V, U = np.meshgrid(v, u)
return np.sqrt(U**2 + V**2)
def aplicar_filtro_freq(img, H):
"""Applies filter H (centered) to image via FFT."""
F = np.fft.fftshift(np.fft.fft2(img.astype(np.float64)))
Fg = F * H
g = np.real(np.fft.ifft2(np.fft.ifftshift(Fg)))
return cv2.normalize(g, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)
M, N = img_gray.shape
D = distancia_centro(M, N)
D0 = 30 # cutoff frequency
n_bw = 2 # Butterworth order
# ── Transfer functions ──────────────────────────────────────────────────
H_ideal = (D <= D0).astype(np.float64)
H_gauss = np.exp(-D**2 / (2 * D0**2))
H_bw = 1.0 / (1.0 + (D / D0)**(2 * n_bw))
# ── Filtered images ─────────────────────────────────────────────────────────
img_ideal = aplicar_filtro_freq(img_gray, H_ideal)
img_gauss = aplicar_filtro_freq(img_gray, H_gauss)
img_bw = aplicar_filtro_freq(img_gray, H_bw)
# ── H(u,v) profiles ─────────────────────────────────────────────────────────
def fig2img(fig):
b = io.BytesIO(); fig.savefig(b, format='png', dpi=100); plt.close(fig); b.seek(0)
return (plt.imread(b)[:,:,:3]*255).astype(np.uint8)
fig, ax = plt.subplots(figsize=(6, 3))
linha = M // 2
ax.plot(H_ideal[linha, :], label="Ideal", color="#D85A30", lw=1.5, ls="--")
ax.plot(H_gauss[linha, :], label="Gaussiano", color="#1D9E75", lw=1.5)
ax.plot(H_bw[linha, :], label="Butterworth", color="#534AB7", lw=1.5)
ax.axvline(N//2-D0, color="#aaa", lw=0.8, ls=":")
ax.axvline(N//2+D0, color="#aaa", lw=0.8, ls=":")
ax.set(title="Perfis H(u,v) — linha central", xlabel="v", ylabel="H(u,v)")
ax.legend(fontsize=8); plt.tight_layout()
perfil_img = fig2img(fig)
# ── Visualized H filters ────────────────────────────────────────────────────
def H_vis(H):
return cv2.normalize((H*255).astype(np.uint8), None, 0, 255, cv2.NORM_MINMAX)
mm.show(
[img_gray, img_ideal, img_gauss, img_bw,
H_vis(H_ideal), H_vis(H_gauss), H_vis(H_bw), perfil_img],
titles=[
"Original", "Ideal LPF", "Gaussian LPF", "Butterworth LPF (n=2)",
"Ideal H", "Gaussian H","Butterworth H", "H(u,v) Profiles"
],
cols=4, figsize=(16, 9)
)
# High-pass filter: complement of the Gaussian low-pass
# Reuses aplicar_filtro_freq() defined in the previous cell
H_alta = 1 - H_gauss
img_alta = aplicar_filtro_freq(img_gray, H_alta)
mm.show(
[img_gray, img_alta],
titles=["Original", "Gaussian high-pass ($D_0=30$)"],
cols=2
)
Periodic noise — associated with electrical interference, regular sensor patterns, or scanning artifacts — appears in the Fourier spectrum as symmetrical point peaks around the center.
The notch filter selectively attenuates these frequencies while preserving the remaining image components. An application example is presented in Figure 5.15.
# ── Image with synthetic periodic noise ──────────────────────────────────────
h_img, w_img = img_gray.shape
x = np.arange(w_img)
y = np.arange(h_img)
X, Y = np.meshgrid(x, y)
# Use integer and consistent frequencies (important for perfect restoration)
u0, v0 = 20, 20 # exact noise frequencies
ruido = 40 * np.sin(2 * np.pi * (u0 * X / w_img + v0 * Y / h_img))
img_ruidosa = np.clip(img_gray.astype(np.float64) + ruido, 0, 255).astype(np.uint8)
# ── Spectrum of the noisy image ──────────────────────────────────────────────
F_r = np.fft.fftshift(np.fft.fft2(img_ruidosa.astype(np.float64)))
mag_r = np.log1p(np.abs(F_r))
mag_vis = cv2.normalize(mag_r, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)
# ── Notch mask ───────────────────────────────────────────────────────────────
mascara = np.ones((h_img, w_img), dtype=np.float64)
r_notch = 8 # notch radius (fine-tune if necessary)
def suprimir_pico(mask, cy, cx, r):
"""Zeros a disk of radius r centered at (cy, cx)"""
yy, xx = np.ogrid[:mask.shape[0], :mask.shape[1]]
dist = np.sqrt((yy - cy)**2 + (xx - cx)**2)
mask[dist <= r] = 0
return mask
# Central coordinates
cy, cx = h_img // 2, w_img // 2
# Suppress the 4 symmetric peaks (important!)
for dy, dx in [(v0, u0), (-v0, -u0), (v0, -u0), (-v0, u0)]:
mascara = suprimir_pico(mascara, cy + dy, cx + dx, r_notch)
mascara_vis = (mascara * 255).astype(np.uint8)
# ── Filtering and reconstruction ─────────────────────────────────────────────
F_filtrada = F_r * mascara
img_rest = np.real(np.fft.ifft2(np.fft.ifftshift(F_filtrada)))
img_rest_vis = cv2.normalize(img_rest, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)
# Evaluation
psnr = cv2.PSNR(img_gray, img_rest_vis)
ssim = cv2.SSIM(img_gray, img_rest_vis) if hasattr(cv2, 'SSIM') else "N/A"
#ssim = ssim_sk(img_gray, img_rest_vis, data_range=255)
print(f"PSNR (original vs restored): {psnr:.2f} dB")
# ── Visualization ────────────────────────────────────────────────────────────
mm.show(
[img_ruidosa, mag_vis, mascara_vis, img_rest_vis],
titles=[
"With periodic noise",
"Spectrum (log)",
"Notch mask",
f"Restored (PSNR={psnr:.1f} dB)"
],
cols=4,
figsize=(16, 4)
)PSNR (original vs restored): 32.93 dB
| Filter | Visual effect | Artifact | Use |
|---|---|---|---|
| Ideal low-pass | Intense smoothing | Ringing | Illustrative |
| Gaussian low-pass | Gentle smoothing | No ringing | General smoothing |
| Butterworth low-pass | Controlled smoothing | Ringing (high orders) | Trade-off between smoothing and selectivity |
| High-pass | Edge enhancement | Noise amplification | Contour detection |
| Notch | Selective frequency removal | Possible local distortions | Periodic noise removal |
The design of filters in the frequency domain consists of defining spectral masks. However, effects in the spatial domain, such as ringing and blurring, emerge directly from these choices in the spectrum.
The Fourier Transform decomposes the signal into global frequencies: each coefficient \(F(u,v)\) receives contributions from the entire image, with no explicit information about the spatial localization of those frequencies. Thus, localized structures, such as edges, are represented in a distributed manner across the spectrum.
Wavelets (ondaletas) overcome this limitation by using basis functions localized in space, which can be translated and scaled. These functions have compact support, that is, they are nonzero only in a finite region of the domain, allowing a simultaneous representation in terms of frequency and spatial localization.
The Fourier Transform accurately describes which frequencies are present in a signal, but it does not explicitly represent where those frequencies occur in space.
In the experiment presented in Figure 5.16, two images with structures located at different positions produce nearly identical magnitude spectra. This occurs because the Fourier representation is global: each coefficient receives contributions from the entire image.
As a result, the magnitude spectrum does not explicitly represent the location of edges or other structures, only the distribution of the frequencies present. This limitation motivated the development of multiresolution representations, such as the Discrete Wavelet Transform (DWT), capable of simultaneously describing the frequency and spatial localization of image structures.
img_sinal1 = np.zeros((128, 128)); img_sinal1[:, 20:25] = 1; img_sinal1[100:105, :] = 1
img_sinal2 = np.zeros((128, 128)); img_sinal2[:, 90:95] = 1; img_sinal2[30:35, :] = 1
mag1 = np.log1p(np.abs(np.fft.fftshift(np.fft.fft2(img_sinal1))))
mag2 = np.log1p(np.abs(np.fft.fftshift(np.fft.fft2(img_sinal2))))
mm.show([img_sinal1, mag1, img_sinal2, mag2],
titles=["Signal A", "Spectrum A", "Signal B (Shifted)", "Spectrum B"],
cols=4, figsize=(14, 4)
)
The Discrete Wavelet Transform (DWT) applies, separately in the horizontal and vertical directions, two complementary filters: a low-pass \(h\) (approximation) and a high-pass \(g\) (details), followed by downsampling by a factor of 2 in each dimension. This process produces four subbands, whose names indicate the combination of filters applied in each direction (L = Low-pass; H = High-pass). The characteristics of each subband are summarized in Table 5.3.
\[ \text{DWT}(f)=\{\underbrace{\text{LL}}_{\text{aprox.}},\; \underbrace{\text{LH}}_{\text{horizontal details}},\; \underbrace{\text{HL}}_{\text{vertical details}},\; \underbrace{\text{HH}}_{\text{diagonal details}}\}. \]
| Subband | Filters applied | Visual content |
|---|---|---|
| LL | low × low | Image approximation (smoothed and reduced version) |
| LH | low × high | Horizontal edges and vertical variations |
| HL | high × low | Vertical edges and horizontal variations |
| HH | high × high | Diagonal details and textures |
The decomposition can be applied recursively to the LL subband, generating a multiresolution representation. After \(J\) levels, a structure with \(3J+1\) subbands is obtained, where each new level reduces the resolution of the approximation component.
The multiresolution decomposition of wavelets has a conceptual relationship with the hierarchical representations used in convolutional neural networks (CNNs). In both cases, successive stages of filtering and resolution reduction produce increasingly abstract descriptions of the image. However, wavelets use mathematically defined and reconstructible filters, whereas CNNs learn their filters during training.
Different wavelet families present distinct trade-offs between spatial support, smoothness, and compression capability. Support corresponds to the extent of the wavelet function in the spatial domain: the smaller the support, the more localized the function; the larger it is, the smoother its representation tends to be, albeit at a higher computational cost. Table 5.4 compares some of the most commonly used families.
| Wavelet | Support length | Vanishing moments | Symmetry | Typical use |
|---|---|---|---|---|
| Haar | 2 | 1 | Asymmetric | Introduction and basic analysis |
| Daubechies db4 | 8 | 4 | Asymmetric | Compression and general analysis |
| Symlet sym4 | 8 | 4 | Nearly symmetric | Signal reconstruction |
| Biorthogonal 5/3 | 5/3 | 2/2 | Symmetric | Lossless JPEG 2000 |
| Biorthogonal 9/7 | 9/7 | 4/4 | Symmetric | Lossy JPEG 2000 |
Vanishing moments measure the wavelet’s ability to represent smooth regions of the image with few nonzero coefficients. A wavelet with \(p\) vanishing moments exactly annihilates polynomials of degree up to \(p-1\). Consequently, the greater the number of vanishing moments, the greater the compression efficiency tends to be in homogeneous regions, although this generally implies functions with longer support.
Figure 5.17 presents the basis functions (wavelets) \(\psi(t)\) in the spatial domain. These functions have compact support, that is, they are nonzero only in a finite region of the domain, unlike the sinusoids of the Fourier Transform, which extend over the entire domain.
# ── Import pywt ───────────────────────────────────────────────────────────────
try:
import pywt
except ImportError:
import subprocess
subprocess.run(["pip", "install", "PyWavelets", "-q"])
import pywt
wavelet_haar = pywt.Wavelet('haar')
wavelet_db4 = pywt.Wavelet('db4')
phi_h, psi_h, x_h = wavelet_haar.wavefun(level=4)
phi_d, psi_d, x_d = wavelet_db4.wavefun(level=4)
fig, ax = plt.subplots(1, 2, figsize=(10, 3))
ax[0].plot(x_h, psi_h, 'b', lw=2); ax[0].set_title("Ondaleta Haar (ψ)")
ax[1].plot(x_d, psi_d, 'g', lw=2); ax[1].set_title("Ondaleta Daubechies 4 (ψ)")
plt.tight_layout(); plt.show()
The diagram in Figure 5.18 illustrates the multiresolution analysis performed by the DWT, in which the approximation subband (LL) is successively decomposed, forming a hierarchical representation with two levels.
The simulator in Figure 5.19 allows for interactive exploration of the 2D Discrete Wavelet Transform (DWT) using the Haar wavelet. The subband decomposition highlights the separation between the approximation component and the detail components of the image.
The different input patterns allow observing the directional behavior of the filters. In images with horizontal and vertical edges, the LH and HL subbands highlight, respectively, the vertical and horizontal intensity variations. In regions of smooth variation, most of the energy is concentrated in the LL approximation subband, while the detail subbands exhibit coefficients close to zero.
Multiresolution analysis can also be observed by increasing the number of decomposition levels. In this case, only the \(\text{LL}_1\) subband is further decomposed, giving rise to the \(\text{LL}_2\), \(\text{LH}_2\), \(\text{HL}_2\), and \(\text{HH}_2\) subbands, which form the second level of the hierarchical representation.
In patterns composed of large homogeneous regions, such as a smooth gradient or a checkerboard made of large blocks, the energy remains predominantly concentrated in the LL subband. In the gradient, this occurs because the differences between neighboring pixels are small. In the checkerboard, in turn, the pixels have practically the same intensity within each block, so that only the boundaries between blocks produce nonzero coefficients in the detail subbands. Since these boundaries occupy only a small fraction of the image, their contribution to the total energy remains reduced.
To enable the visual analysis of these subtle variations, the simulator incorporates a contrast gain control for the details (ranging from 1 to 8). This parameter functions as a linear amplification factor applied exclusively to the coefficients of the detail subbands (LH, HL, and HH) before their on-screen rendering. In scenarios of smooth transition (such as the gradient) or local uniformity (such as the interior of the checkerboard blocks), the numerical differences calculated by the Haar high-pass filter result in coefficients very close to zero, which would render the corresponding quadrants dark and imperceptible to the naked eye. By multiplying these values by the gain, the simulator visually recovers the hidden high-frequency structures and enhances the orientation of the remaining edges.
The energy-per-subband plot quantifies this distribution between the approximation component and the detail components, demonstrating that the visual gain does not alter the original energy metric. In natural images, most of the energy is concentrated in the LL subband, while the LH, HL, and HH subbands mainly represent edges, textures, and other local intensity variations.
The 2D Discrete Wavelet Transform (DWT) decomposes an image into approximation and detail components, organized hierarchically across different scales and orientations. As detail subbands in natural images often exhibit low-contrast coefficients, the practical examples below use a synthetic geometric pattern generated in Python. This approach replicates the behavior of the simulator in Figure 5.19, making the effects of spatial filtering and multiresolution decomposition visually explicit.
Figure 5.20 illustrates the hierarchical structure of the two-level DWT using the Haar wavelet. The process is based on the combined application of low-pass and high-pass filters in the horizontal and vertical directions, followed by subsampling by a factor of 2.
In the first level, the original image yields the approximation subband (\(LL_1\)) and the horizontal (\(LH_1\)), vertical (\(HL_1\)), and diagonal (\(HH_1\)) detail components. In multiresolution analysis, the \(LL_1\) subband is again filtered and subsampled, generating the second decomposition level (\(LL_2\), \(LH_2\), \(HL_2\), and \(HH_2\)).
To enable the visual interpretation of the detail components, the code extracts the absolute value of their coefficients and applies linear (min-max) normalization to occupy the entire dynamic range of gray levels [0, 255]. This operation transforms homogeneous regions (zero coefficients) into black and highlights in white the edges and textures extracted at each scale and orientation.
try:
import pywt
HAS_PYWT = True
except ImportError:
import subprocess
subprocess.run(["pip", "install", "PyWavelets", "-q"])
import pywt
HAS_PYWT = True
import numpy as np
import cv2
# ── Generation of the Synthetic Image (Same 'combined' pattern as the simulator) ────────
def gerar_imagem_sintetica(N=256):
img = np.zeros((N, N), dtype=np.float64)
for y in range(N):
for x in range(N):
v = 55 + 35 * (x / N) + 15 * np.sin(y / 24)
# Square
if 24 < x < 100 and 24 < y < 100:
v = 225
# Circle
cx, cy, r = 190, 76, 34
if (x - cx)**2 + (y - cy)**2 < r**2:
v = 205
# Periodic texture (bottom)
if y > 164 and y < 244:
p = 12
v = 185 if ((x // p + y // p) % 2 == 0) else 65
# Diagonal line
if abs(x - y) < 4:
v = 240
img[y, x] = np.clip(v, 0, 255)
return img.astype(np.uint8)
# Replaces the dark coins image with the clear synthetic pattern
img_gray = gerar_imagem_sintetica(256)
# ── 2-level wavelet decomposition ─────────────────────────────────────────────
wavelet = "haar"
img_float = img_gray.astype(np.float64)
# Level 1
coefs1 = pywt.dwt2(img_float, wavelet)
LL1, (LH1, HL1, HH1) = coefs1
# Level 2 (applied on LL1)
coefs2 = pywt.dwt2(LL1, wavelet)
LL2, (LH2, HL2, HH2) = coefs2
def sb_vis(sb):
"""Normalizes subband for visualization [0,255]."""
return cv2.normalize(np.abs(sb), None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)
print(f"Original shape : {img_gray.shape}")
print(f"LL1 (level 1) : {LL1.shape} | LH1/HL1/HH1: {LH1.shape}")
print(f"LL2 (level 2) : {LL2.shape} | LH2/HL2/HH2: {LH2.shape}")
imgs_dwt = [img_gray, sb_vis(LL1), sb_vis(LH1), sb_vis(HL1), sb_vis(HH1),
sb_vis(LL2), sb_vis(LH2), sb_vis(HL2), sb_vis(HH2)]
titles_dwt = ["Original",
"LL₁ (aprox.)", "LH₁ (horiz.)", "HL₁ (vert.)", "HH₁ (diag.)",
"LL₂ (aprox.)", "LH₂ (horiz.)", "HL₂ (vert.)", "HH₂ (diag.)"]
mm.show(imgs_dwt, titles=titles_dwt, cols=5, figsize=(16, 7))Original shape : (256, 256)
LL1 (level 1) : (128, 128) | LH1/HL1/HH1: (128, 128)
LL2 (level 2) : (64, 64) | LH2/HL2/HH2: (64, 64)
The choice of the basis function (wavelet) directly influences how image features are distributed and encoded by the DWT coefficients. Figure 5.21 compares the practical results obtained by applying four distinct families to the synthetic geometric pattern: haar, db4, sym4, and bior2.2.
Due to its short support and step-function shape, the Haar wavelet produces coefficients that are highly localized at spatial discontinuities, generating thin, sharp edges in the detail subbands. In contrast, families such as Daubechies (db4) and Symlets (sym4), which exhibit larger support (longer filters) and a greater number of vanishing moments, yield smoother, more distributed responses around transitions, which can introduce slight oscillations or blurring at abrupt boundaries.
This behavior highlights the classic trade-off in multiresolution analysis: smaller supports favor precise spatial localization of edges, while larger supports and a greater number of vanishing moments tend to produce sparser, smoother representations. This smoothness and ability to attenuate high frequencies ensure greater efficiency in energy compaction—fundamental characteristics for data compression and denoising applications.
import numpy as np
import cv2
import pywt
# Ensures that img_gray and img_float use the same clear synthetic pattern
if 'gerar_imagem_sintetica' in globals():
img_gray = gerar_imagem_sintetica(256)
else:
# Fallback in case the previous block was not executed in the same session
def gerar_imagem_sintetica(N=256):
img = np.zeros((N, N), dtype=np.float64)
for y in range(N):
for x in range(N):
v = 55 + 35 * (x / N) + 15 * np.sin(y / 24)
if 24 < x < 100 and 24 < y < 100: v = 225
cx, cy, r = 190, 76, 34
if (x - cx)**2 + (y - cy)**2 < r**2: v = 205
if y > 164 and y < 244:
p = 12
v = 185 if ((x // p + y // p) % 2 == 0) else 65
if abs(x - y) < 4: v = 240
img[y, x] = np.clip(v, 0, 255)
return img.astype(np.uint8)
img_gray = gerar_imagem_sintetica(256)
img_float = img_gray.astype(np.float64)
wavelets_comp = ["haar", "db4", "sym4", "bior2.2"]
imgs_comp, titles_comp = [], []
for wname in wavelets_comp:
LL, (LH, HL, HH) = pywt.dwt2(img_float, wname)
imgs_comp += [sb_vis(LL), sb_vis(HH)]
titles_comp += [f"{wname} — LL₁", f"{wname} — HH₁"]
mm.show(imgs_comp, titles=titles_comp, cols=4, figsize=(14, 8))
One of the main applications of the Discrete Wavelet Transform (DWT) is data compression, driven by the sparse representation capability of the coefficients. Figure 5.22 illustrates the effect of hard thresholding, a technique in which detail coefficients with magnitude below a threshold \(T\) are entirely zeroed out before the synthesis process carried out by the Inverse Discrete Wavelet Transform (IDWT).
As the threshold \(T\) is increased, a growing number of high-frequency coefficients are zeroed. Since they concentrate less energy, removing these components considerably reduces the amount of information required to represent the image, while keeping the global approximation component (the deepest \(LL\) subband) intact to preserve the macro structure. Visually, this discarding of coefficients manifests itself through the progressive disappearance of fine textures and the smoothing of abrupt intensity transitions.
The fidelity of the reconstructed image relative to the original is quantified by the Peak Signal-to-Noise Ratio (PSNR) metric, expressed in decibels (dB). Higher PSNR values indicate less distortion and greater mathematical proximity to the original signal. The practical experiment highlights the gradual decay of PSNR as the thresholding aggressiveness increases, allowing the optimal threshold for the balance between compression and visual degradation to be numerically evaluated.
import numpy as np
import cv2
import pywt
# Ensures img_gray uses the same clear synthetic pattern
if 'gerar_imagem_sintetica' in globals():
img_gray = gerar_imagem_sintetica(256)
else:
# Fallback in case the previous block was not executed in the same session
def gerar_imagem_sintetica(N=256):
img = np.zeros((N, N), dtype=np.float64)
for y in range(N):
for x in range(N):
v = 55 + 35 * (x / N) + 15 * np.sin(y / 24)
if 24 < x < 100 and 24 < y < 100: v = 225
cx, cy, r = 190, 76, 34
if (x - cx)**2 + (y - cy)**2 < r**2: v = 205
if y > 164 and y < 244:
p = 12
v = 185 if ((x // p + y // p) % 2 == 0) else 65
if abs(x - y) < 4: v = 240
img[y, x] = np.clip(v, 0, 255)
return img.astype(np.uint8)
img_gray = gerar_imagem_sintetica(256)
def dwt_threshold_reconstruct(img, wavelet='db4', nivel=2, threshold=0.0):
"""Decomposes, applies threshold and reconstructs via IDWT."""
coefs = pywt.wavedec2(img.astype(np.float64), wavelet, level=nivel)
# Copies and applies hard thresholding to all details
coefs_t = [coefs[0]] # Final LL is not thresholded
for detalhe in coefs[1:]:
coefs_t.append(tuple(pywt.threshold(sb, threshold, mode='hard') for sb in detalhe))
rec = pywt.waverec2(coefs_t, wavelet)
# Crop to original dimension
rec = rec[:img.shape[0], :img.shape[1]]
return np.clip(rec, 0, 255).astype(np.uint8)
thresholds = [0, 10, 30, 60, 100]
imgs_thr = [img_gray]
titles_thr = ["Original"]
for t in thresholds:
rec = dwt_threshold_reconstruct(img_gray, threshold=t)
psnr = cv2.PSNR(img_gray, rec)
imgs_thr.append(rec)
titles_thr.append(f"T={t} PSNR={psnr:.1f} dB")
mm.show(imgs_thr, titles=titles_thr, cols=3, figsize=(14, 10))
Table 5.5 summarizes the main structural and operational differences between the Discrete Fourier Transform (DFT) and the Discrete Wavelet Transform (DWT).
| Criterion | Fourier (DFT) | Wavelet (DWT) |
|---|---|---|
| Basis functions | Sinusoids with infinite support | Functions with compact support |
| Spatial localization | Not explicit (global) | Explicit (local) |
| Spectral filtering | Excellent for fine frequency control | Subband-based (scales) |
| Image compression | Basis of DCT (traditional JPEG) | Basis of DWT (JPEG 2000) |
| Multiscale analysis | No | Yes |
| Periodic noise removal | Highly efficient | Not recommended |
| Non-stationary signals | Limited | Highly efficient |
In practical terms, the DFT stands out as the ideal tool for pure spectral analysis, design of selective filters in the frequency domain, and attenuation of periodic and harmonic noises. On the other hand, the DWT excels in scenarios that require rigorous preservation of the spatial localization of features associated with their frequency content, standing out in data compression, multiresolution analysis, and processing of abrupt transitions. Thus, both transforms should be understood as perfectly complementary techniques, mapping distinct and specific paths for solving problems in DIP-CV.
When drawing analogies between image and audio processing, it is important to consider the fundamental differences:
In stereo/multichannel audio systems, the phase between channels is crucial for the perception of spatial location (interaural phase and time differences).
In monaural systems, phase has limited perceptual influence — the human ear is relatively insensitive to the absolute phase of isolated sinusoidal components.
In images, the DFT phase is always fundamental for the spatial location of structures, regardless of whether it is a monochromatic or color image.
The analogy between phase in audio and phase in images should be used with caution, emphasizing that, although both carry information about the spatial/temporal organization of the signal, the perceptual mechanisms are fundamentally different.
While wavelets establish the theoretical foundation of the JPEG 2000 standard, the traditional JPEG standard is based on the Discrete Cosine Transform (DCT). Despite structural differences, both approaches share the same fundamental principle: compacting image energy into a reduced number of coefficients and discarding the least relevant components with minimal visual impact.
The central goal of compression is to reduce the volume of data required for storing or transmitting an image. This process is made feasible by identifying and eliminating structural and perceptual redundancies.
The development of compression algorithms is based on the identification and elimination of three main categories of redundancy, summarized in Table 5.6.
| Type | Definition | Exploration Approach |
|---|---|---|
| Spatial (interpixel) | High correlation and statistical dependence between neighboring pixels. | DCT, DWT, and predictive coding. |
| Spectral (interchannel) | Statistical correlation between the color channels of the same image. | Color space transformations (e.g., RGB to \(YC_bC_r\)). |
| Psychovisual | Insensitivity of the human visual system (HVS) to high-frequency and low-contrast variations. | Selective coefficient quantization processes. |
Depending on the preservation of the original information after the decoding process, compression methods are divided into two fundamental classes:
The Discrete Cosine Transform (DCT) constitutes the central operation of the JPEG standard. Unlike the DFT, which uses a complex basis, the DCT is based on purely real trigonometric functions. For an image block \(f(x,y)\) of dimensions \(N \times N\), the 2D DCT-II maps the spatial signal to the spatial frequency domain, generating the coefficient matrix \(C(u,v)\) through:
\[ C(u,v) = \alpha(u)\,\alpha(v) \sum_{x=0}^{N-1}\sum_{y=0}^{N-1} f(x,y)\, \cos\!\left[\frac{\pi(2x+1)u}{2N}\right] \cos\!\left[\frac{\pi(2y+1)v}{2N}\right] \tag{5.9}\]
where the orthogonal normalization factors are given by \(\alpha(0) = \sqrt{1/N}\) and \(\alpha(k) = \sqrt{2/N}\) for \(k > 0\).
Each coefficient \(C(u,v)\) quantifies the contribution — or “weight” — of a specific spatial frequency within that block. The term \(C(0,0)\) is called the DC component and represents the average intensity of the block (zero frequency). The remaining coefficients, known as AC components (Alternating Current), correspond to progressively higher spatial frequencies.
From a geometric perspective, Equation 5.9 performs the projection of the pixel block onto a set of orthogonal functions. For the standard JPEG case (\(N=8\)), the spatial block is decomposed into a linear combination of 64 two-dimensional basis functions, denoted by \(B_{u,v}(x,y)\) and generated by the product of cosine functions:
\[B_{u,v}(x,y) = \cos\left[ \frac{\pi (2x+1)u}{16} \right] \cos\left[ \frac{\pi (2y+1)v}{16} \right]\]
Thus, the inverse operation can be interpreted as the exact reconstruction of the original block through the weighted sum of these 64 basis matrices, where each coefficient \(C(u,v)\) acts as the analytical weight of its respective harmonic component.
The spatial frequency indicated by the indices \((u,v)\) determines the number of oscillation cycles along the horizontal and vertical dimensions of the block. As illustrated in Figure 5.23 — whose code isolates each basis by applying the inverse transformation to unit impulses —, these 64 functions are organized into an \(8 \times 8\) matrix. The upper-left corner (\(u=0, v=0\)) displays the uniform zero-frequency (DC) pattern, while moving to the right (the \(u\) axis) or downward (the \(v\) axis) maps progressively larger harmonic variations, representing rapid transitions, edges, and textures in the horizontal, vertical, and diagonal orientations.
Both the DCT and the DFT map a spatial \(N \times N\) block into a coefficient matrix of the same dimensions. However, for natural images, the DCT exhibits greater efficiency in energy compaction at low frequencies. This occurs because the DCT implicitly assumes an even symmetry of the signal at the block boundaries, which is equivalent to a continuous periodic extension, minimizing the effect of spectral spreading (ringing). As a result, most AC coefficients decay rapidly to values close to zero, optimizing the compression pipeline without introducing perceptible visual degradation.
from scipy.fft import dct, idct # added line
fig, axes = plt.subplots(8, 8, figsize=(6, 6))
fig.subplots_adjust(hspace=0.05, wspace=0.05)
for i in range(8):
for j in range(8):
coef = np.zeros((8, 8)); coef[i, j] = 1
b = idct(idct(coef.T, norm='ortho').T, norm='ortho')
axes[i, j].imshow(b, cmap='gray')
axes[i, j].axis('off')
plt.suptitle("As 64 Bases da DCT 8x8", y=0.92, fontsize=12, fontweight='bold')
plt.show()
Before applying the DCT, the pixels of the intensity block are routinely shifted (by subtracting \(128\) for 8-bit images) in order to center the signal around zero, eliminating unnecessary DC components. When computing the DCT on the resulting block, the property of energy compaction becomes evident: almost all of the variance and information from the original image is concentrated in the DC coefficient (\(C(0,0)\)) and in the first low-frequency AC harmonics.
Figure 5.24 demonstrates this phenomenon through a progressive reconstruction by abrupt truncation. Instead of using all 64 coefficients, the algorithm preserves only the first \(k\) components—selected based on a scan that prioritizes low spatial frequencies—and nullifies the remaining ones.
The inverse synthesis (IDCT) performed with only a fraction of the coefficients (such as 15% or 30%) is already capable of recovering the structures and the macro illumination of the original pixel block. As higher-frequency harmonics are progressively reintroduced, fine details and rapid transitions are restored. This behavior validates the principle of perceptual compression: the discarded high frequencies carry little energy, and their absence, under normal conditions, has a secondary visual impact on the observer’s perception.
from scipy.fft import dct, idct
def dct2(bloco):
"""Orthogonal 2D DCT-II (separable)."""
return dct(dct(bloco.T, norm='ortho').T, norm='ortho')
def idct2(coefs):
"""Orthogonal 2D IDCT-II."""
return idct(idct(coefs.T, norm='ortho').T, norm='ortho')
# ── Centralized 8×8 block of the image ─────────────────────────────────────────
cy, cx = img_gray.shape[0]//2, img_gray.shape[1]//2
bloco = img_gray[cy:cy+8, cx:cx+8].astype(np.float64) - 128.0
C = dct2(bloco)
print("DCT coefficients of the 8×8 block:")
print(np.round(C).astype(int))
print(f"\nDC energy : {C[0,0]**2:.1f}")
print(f"Total energy : {(C**2).sum():.1f}")
print(f"DC fraction : {C[0,0]**2 / (C**2).sum():.1%} ← energy concentration")
# ── Progressive reconstruction ──────────────────────────────────────────────────
imgs_rec = [cv2.normalize((bloco+128).astype(np.uint8), None, 0, 255, cv2.NORM_MINMAX)]
titles_rec = ["Bloco original\n(8×8 pixels)"]
for keep in [1, 4, 10, 20, 40, 64]:
C_trunc = np.zeros_like(C)
indices = sorted([(u,v) for u in range(8) for v in range(8)], key=lambda p: p[0]+p[1])
for u, v in indices[:keep]:
C_trunc[u, v] = C[u, v]
rec = np.clip(idct2(C_trunc) + 128, 0, 255).astype(np.uint8)
imgs_rec.append(rec)
titles_rec.append(f"{keep} coef.\n({keep/64:.0%} do total)")
mm.show(imgs_rec, titles=titles_rec, cols=4, figsize=(12, 7))DCT coefficients of the 8×8 block:
[[ 450 2 -199 -1 0 0 -14 0]
[ -2 506 0 -58 0 24 0 -24]
[-199 2 89 -1 82 1 0 0]
[ -1 -58 0 -158 0 59 0 -24]
[ 0 0 83 0 -89 0 -34 0]
[ 0 24 0 58 0 84 0 -58]
[ -14 0 0 0 -34 0 89 -1]
[ 0 -24 0 -24 0 -58 0 -76]]
DC energy : 202725.1
Total energy : 639868.0
DC fraction : 31.7% ← energy concentration
The JPEG standard operates by dividing the image into disjoint blocks of \(8 \times 8\) pixels, processed through a sequence of spatial, perceptual, and statistical transformations. The complete encoding pipeline is structured into six main stages:
\[ \text{RGB} \xrightarrow{\text{(1) } YC_bC_r} \xrightarrow{\text{(2) Subsampling}} \xrightarrow{\text{(3) Blocks } 8 \times 8} \xrightarrow{\text{(4) DCT}} \xrightarrow{\text{(5) Quantization}} \xrightarrow{\text{(6) Entropy Coding}} \]
Table 5.7 details the analytical function and the perceptual rationale that justifies each of these stages.
| Stage | Operation | Perceptual and Statistical Rationale |
|---|---|---|
| 1 | Conversion \(RGB \rightarrow YC_bC_r\) | Separates luminance (\(Y\)) from chrominance (\(C_b, C_r\)). The human visual system (HVS) exhibits greater sensitivity to variations in brightness than to color. |
| 2 | Chrominance subsampling (e.g., 4:2:0) | Reduces the spatial resolution of the color channels by half, discarding redundant data with negligible visual impact. |
| 3–4 | Centering and application of the \(8 \times 8\) DCT | Translates the pixels to the range \([-128, 127]\) and compacts the spectral energy of the block into the low-frequency coefficients. |
| 5 | Selective linear quantization | Divides each coefficient \(C(u,v)\) by the corresponding element of the matrix \(Q(u,v)\), applying integer rounding. This constitutes the primary source of lossy compression. |
| 6 | Zigzag scanning and coding | Orders the quantized coefficients to maximize consecutive null sequences, optimizing run-length encoding (RLE) and Huffman coding. |
The quantization matrix \(Q(u,v)\) is the central mechanism for controlling the trade-off between compression ratio and visual quality. In the practical algorithm of Figure 5.25, the quality factor specified by the user (on a scale from 1 to 100) is converted into a scalar that parameterizes the severity of the matrix \(Q\). Reduced quality values expand the divisors of \(Q(u,v)\), forcing the massive truncation of AC coefficients to zero. When this elimination is excessive, the discontinuity at the boundaries of adjacent blocks is not attenuated during reconstruction, giving rise to the so-called blocking artifacts.
The efficiency of the entropy coder subsequent to quantization depends directly on the ordering of the data. Since the DCT concentrates the essential energy at the upper-left vertex of the matrix (low frequencies) and pushes the null coefficients toward the opposite ends, a linear reading by rows or columns would fragment the sequences of zeros.
Zigzag ordering solves this limitation by traversing the matrix diagonally in increasing order of spatial frequency. This mapping groups the significant coefficients at the beginning of the vector and concentrates the null coefficients into a single continuous sequence at the end of the arrangement, allowing the RLE algorithm to encode large blocks of data compactly and efficiently.
RLE (Run-Length Encoding) is a lossless compression technique that encodes consecutive sequences of identical values — especially zeros — as a pair (count, value). In JPEG, after zigzag scanning, the quantized coefficients are organized so that the zeros are concentrated at the end of the vector. RLE then compresses this long run of zeros with extreme efficiency, optimizing the storage and transmission of the compressed image.
import numpy as np
import cv2
from scipy.fft import dct, idct
# ── Safe Loading of the Camera Image (skimage) ────────────────────────────────
try:
from skimage import data
img_gray = data.camera()
except ImportError:
import subprocess
subprocess.run(["pip", "install", "scikit-image", "-q"])
from skimage import data
img_gray = data.camera()
# Slightly resizes to 256x256 to keep the pattern and speed of the previous tests
img_gray = cv2.resize(img_gray, (256, 256))
# ── Luminance quantization table (standard JPEG) ──────────────────────────────
Q_luma = np.array([
[16,11,10,16,24,40,51,61],
[12,12,14,19,26,58,60,55],
[14,13,16,24,40,57,69,56],
[14,17,22,29,51,87,80,62],
[18,22,37,56,68,109,103,77],
[24,35,55,64,81,104,113,92],
[49,64,78,87,103,121,120,101],
[72,92,95,98,112,100,103,99]
], dtype=np.float64)
def dct2(bloco):
"""2D orthogonal DCT-II (separable)."""
return dct(dct(bloco.T, norm='ortho').T, norm='ortho')
def idct2(coefs):
"""2D orthogonal IDCT-II."""
return idct(idct(coefs.T, norm='ortho').T, norm='ortho')
def jpeg_compress_block(bloco, Q_table):
"""DCT → quantization → dequantization → IDCT in 8×8 block."""
C = dct2(bloco.astype(np.float64) - 128)
Cq = np.round(C / Q_table) * Q_table # quantizes and dequantizes
return np.clip(idct2(Cq) + 128, 0, 255)
def jpeg_quality_compress(img, qualidade=50):
"""Simplified JPEG: compresses entire image by 8×8 blocks."""
if qualidade < 50:
escala = 5000 / qualidade
else:
escala = 200 - 2 * qualidade
# Fixed from 'scala' to 'escala'
Q = np.clip(np.round(Q_luma * escala / 100), 1, 255)
h, w = img.shape
result = np.zeros_like(img, dtype=np.float64)
for r in range(0, h-7, 8):
for c in range(0, w-7, 8):
result[r:r+8, c:c+8] = jpeg_compress_block(img[r:r+8, c:c+8], Q)
return result.astype(np.uint8)
# ── Comparison of quality factors ─────────────────────────────────────────────
qualidades = [10, 25, 50, 75, 90]
imgs_jpeg = [img_gray]
titles_jpeg = ["Original\n(Cameraman)"]
for q in qualidades:
rec = jpeg_quality_compress(img_gray, qualidade=q)
psnr = cv2.PSNR(img_gray, rec)
imgs_jpeg.append(rec)
titles_jpeg.append(f"Q={q}\nPSNR={psnr:.1f}dB")
mm.show(imgs_jpeg, titles=titles_jpeg, cols=3, figsize=(14, 10))
The simulator in Figure 5.26 allows exploring the impact of the quantization process on an \(8 \times 8\) block extracted from a real image, synthesizing in real time the following components:
The choice of a digital storage format directly impacts the trade-off between visual quality, file size, and computational decoding cost. The three most relevant formats for web architectures and visual computing systems are JPEG, PNG, and WebP.
Table 5.8 synthesizes the structural properties of the main rasterized image formats.
| Characteristic | JPEG | PNG | WebP |
|---|---|---|---|
| Compression | Lossy | Lossless | Lossy and lossless. |
| Transparency (alpha channel) | No | Yes | Yes. |
| Animation support | No | Limited (APNG) | Yes. |
| Base algorithm | DCT + Huffman | DEFLATE (LZ77 + Huffman) | VP8 / VP8L. |
| Best for | Photography | Graphics, text, and icons | Universal use in a Web environment. |
| Worst for | Text and sharp edges | Complex photographic images | Legacy compatibility. |
Two objective metrics are widely adopted to quantify the distortion introduced by compression processes:
Peak Signal-to-Noise Ratio (PSNR): \[ \text{PSNR} = 10\,\log_{10}\!\left(\frac{L^2}{\text{MSE}}\right) \quad [\text{dB}] \tag{5.10}\]
where \(L = 255\) for 8-bit quantized images and \(\text{MSE}\) represents the Mean Squared Error. PSNR values above 40 dB indicate excellent fidelity; between 30 dB and 40 dB represent good quality; and values below 30 dB correspond to easily perceptible visual degradations.
Structural Similarity Index (SSIM): \[ \text{SSIM}(f,g) = \frac{(2\mu_f\mu_g + c_1)(2\sigma_{fg} + c_2)}{(\mu_f^2+\mu_g^2+c_1)(\sigma_f^2+\sigma_g^2+c_2)} \tag{5.11}\]
SSIM evaluates local windows of the image based on three complementary components: luminance (\(\mu_f, \mu_g\)), contrast (\(\sigma_f, \sigma_g\)), and structure (\(\sigma_{fg}\)), weighted by stability constants \(c_1\) and \(c_2\). The index ranges over the interval \([-1, 1]\), where unity represents perfect identity. Unlike PSNR, SSIM considers the spatial organization of errors, aligning with the perception of the human visual system (HVS).
PSNR has a simple mathematical formulation and low computational cost; however, it tends to overestimate quality in images with localized distortions or underestimate it in global brightness variations tolerated by the observer. SSIM models biological perception with greater fidelity but requires greater processing effort. For rigorous analyses of codecs, it is recommended to report both statistical metrics in a complementary manner.
The mathematical nature of the codec dictates the type of degradation introduced at reduced bit rates. As illustrated in Figure 5.27, aggressive DCT-based compression in the JPEG standard segments the image into rigid grids, generating blocking artifacts. In contrast, algorithms based on predictive coding or representations subjected to advanced spatial transforms (such as WebP and JPEG 2000) eliminate block discontinuities but introduce loss of fine texture and characteristic blurring around high-contrast edges.
import os
import cv2
# Ensures the directory exists and saves the compressed files
os.makedirs("imagens/comp_test", exist_ok=True)
cv2.imwrite("imagens/comp_test/camera_q10.jpg", img_gray, [cv2.IMWRITE_JPEG_QUALITY, 10])
cv2.imwrite("imagens/comp_test/camera_q10.webp", img_gray, [cv2.IMWRITE_WEBP_QUALITY, 10])
# Extraction of region of interest for artifact visualization (4x zoom)
zoom_original = cv2.resize(img_gray[120:200, 150:230], (320, 320),
interpolation=cv2.INTER_NEAREST)
rec_jpeg = cv2.imread("imagens/comp_test/camera_q10.jpg", cv2.IMREAD_GRAYSCALE)
zoom_jpeg = cv2.resize(rec_jpeg[120:200, 150:230], (320, 320),
interpolation=cv2.INTER_NEAREST)
rec_webp = cv2.imread("imagens/comp_test/camera_q10.webp", cv2.IMREAD_GRAYSCALE)
zoom_webp = cv2.resize(rec_webp[120:200, 150:230], (320, 320),
interpolation=cv2.INTER_NEAREST)
mm.show([zoom_original, zoom_jpeg, zoom_webp],
titles=["Original Zoom", "JPEG Q=10 (Block Artifact)", "WebP Q=10 (Smoothing)"],
cols=3, figsize=(14, 5))
Validation of lossy compression algorithms requires an analysis that correlates storage cost with the fidelity of the reconstructed signal. This assessment is carried out complementarily through global performance curves and the local mapping of distortions induced by the coders.
Figure 5.28 presents the empirical evaluation of the JPEG and WebP pipeline through rate-distortion curves, which monitor compression gain (file size in KB) as a function of PSNR. The PNG format serves as an ideal baseline (\(\text{PSNR} = \infty\)), as its lossless nature prevents any degradation, although it demands a substantially larger data volume.
Analysis of the curves demonstrates the superiority and efficiency of the WebP standard over traditional JPEG: to achieve the same level of mathematical fidelity (such as the excellent quality range, where \(\text{PSNR} > 40\text{ dB}\)), the WebP encoder generates significantly smaller files. This behavior reflects the practical impact of algorithmic evolution on optimizing digital transmission and storage systems.
import os
import cv2
import matplotlib.pyplot as plt
# Ensures the test directory exists
os.makedirs("imagens/comp_test", exist_ok=True)
resultados = []
# ── JPEG ──────────────────────────────────────────────────────────────────────
for q in [10, 20, 30, 40, 50, 60, 70, 80, 90, 95]:
path = f"imagens/comp_test/camera_q{q}.jpg"
cv2.imwrite(path, img_gray, [cv2.IMWRITE_JPEG_QUALITY, q])
rec = cv2.imread(path, cv2.IMREAD_GRAYSCALE)
resultados.append({"formato": "JPEG", "qualidade": q,
"PSNR": cv2.PSNR(img_gray, rec),
"KB": os.path.getsize(path)/1024})
# ── PNG ───────────────────────────────────────────────────────────────────────
path_png = "imagens/comp_test/camera.png"
cv2.imwrite(path_png, img_gray, [cv2.IMWRITE_PNG_COMPRESSION, 9])
resultados.append({"formato": "PNG", "qualidade": "lossless",
"PSNR": float('inf'), "KB": os.path.getsize(path_png)/1024})
# ── WebP ──────────────────────────────────────────────────────────────────────
for q in [50, 75, 90]:
path_w = f"imagens/comp_test/camera_q{q}.webp"
cv2.imwrite(path_w, img_gray, [cv2.IMWRITE_WEBP_QUALITY, q])
rec_w = cv2.imread(path_w, cv2.IMREAD_GRAYSCALE)
resultados.append(
{"formato": "WebP", "qualidade": q,
"PSNR": cv2.PSUB_VAL if 'cv2.PSNR' in globals() else cv2.PSNR(img_gray, rec_w),
"KB": os.path.getsize(path_w)/1024})
# ── Rate-Distortion Curve Generation ──────────────────────────────────────────
jpeg_r = [r for r in resultados if r["formato"]=="JPEG"]
webp_r = [r for r in resultados if r["formato"]=="WebP"]
png_r = [r for r in resultados if r["formato"]=="PNG"]
fig, ax = plt.subplots(figsize=(8, 4.5))
ax.plot([r["KB"] for r in jpeg_r], [r["PSNR"] for r in jpeg_r],
"o-", label="JPEG", color="#D85A30", lw=2, ms=5)
ax.plot([r["KB"] for r in webp_r], [r["PSNR"] for r in webp_r],
"s-", label="WebP", color="#534AB7", lw=2, ms=5)
ax.axhline(50, color="#1D9E75", lw=2, ls="--",
label=f"PNG sem perda ({png_r[0]['KB']:.1f} KB)")
ax.axhspan(40, 60, alpha=0.05, color="#1D9E75", label="Qualidade excelente (PSNR>40)")
ax.set(xlabel="Tamanho do arquivo (KB)", ylabel="PSNR (dB)",
title="Curva Taxa-Distorção: JPEG × WebP × PNG")
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"\nRaw size (uncompressed): {img_gray.nbytes/1024:.0f} KB")
print(f"\n{'Format':>8} {'Qual.':>6} {'KB':>7} {'PSNR (dB)':>11}")
print("-"*38)
for r in resultados:
psnr_s = f"{r['PSNR']:>11.2f}" if r['PSNR']!=float('inf') else f"{'∞ (lossless)':>11}"
print(f"{r['formato']:>8} {str(r['qualidade']):>6} {r['KB']:>7.1f} {psnr_s}")
Raw size (uncompressed): 64 KB
Format Qual. KB PSNR (dB)
--------------------------------------
JPEG 10 2.3 28.00
JPEG 20 3.5 30.09
JPEG 30 4.5 31.23
JPEG 40 5.4 32.10
JPEG 50 6.2 32.81
JPEG 60 7.0 33.49
JPEG 70 8.5 34.48
JPEG 80 10.7 36.19
JPEG 90 15.7 40.02
JPEG 95 22.3 44.77
PNG lossless 36.2 ∞ (lossless)
WebP 50 4.7 34.29
WebP 75 6.7 36.69
WebP 90 12.5 43.14The Cameraman image (\(256 \times 256\) pixels in grayscale) occupies 64 KB in raw format (uncompressed). As a reference, lossless PNG compresses this volume to 36.2 KB — highlighting that lossless compression already significantly reduces storage for images with homogeneous regions. In contrast, lossy formats (JPEG and WebP) achieve even smaller sizes: JPEG at quality 95 occupies 22.3 KB (PSNR ≈ 45 dB), while WebP at quality 90 reaches 12.5 KB with equivalent PSNR, demonstrating its superior compression efficiency.
Although PSNR provides a rapid numerical indicator, global metrics fail to discriminate how information loss is geometrically distributed across the image. Figure 5.29 addresses this limitation by associating reconstructions at different qualities with their respective absolute error and SSIM maps.
The residual maps—obtained by the normalized absolute difference between the original and compressed images—reveal the intrinsic spatial signature of each coding architecture:
SSIM captures this morphological degradation far more sensitively than PSNR, penalizing the final score as structural organization and fine textures—to which the human visual system is highly responsive—are eliminated by the encoder.
import os
import numpy as np
import cv2
try:
from skimage.metrics import structural_similarity as ssim
except ImportError:
import subprocess
subprocess.run(["pip", "install", "scikit-image", "-q"])
from skimage.metrics import structural_similarity as ssim
# Ensures the tests directory exists
os.makedirs("imagens/comp_test", exist_ok=True)
qualidades_ssim = [25, 50, 75, 95]
imgs_ssim = [img_gray]
titles_ssim = ["Original"]
for q in qualidades_ssim:
path = f"imagens/comp_test/camera_ssim_q{q}.jpg"
# FORCED RECORDING: Generates and records the JPEG with current quality at the correct path
img_compactada = jpeg_quality_compress(img_gray, qualidade=q)
cv2.imwrite(path, img_compactada)
# Safe reading of the newly recorded file
rec = cv2.imread(path, cv2.IMREAD_GRAYSCALE)
if rec is None:
continue
if rec.shape != img_gray.shape:
rec = cv2.resize(rec, (img_gray.shape[1], img_gray.shape[0]))
psnr_v = cv2.PSNR(img_gray, rec)
ssim_v, _ = ssim(img_gray, rec, full=True)
# Normalized absolute difference to highlight the spatial structure of the error
diff_vis = cv2.normalize(np.abs(img_gray.astype(float) - rec.astype(float)),
None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)
imgs_ssim += [rec, diff_vis]
titles_ssim += [f"Q={q}\nPSNR={psnr_v:.1f}dB | SSIM={ssim_v:.3f}",
f"Mapa de erro (Q={q})\n(Bordas e blocagem)"]
mm.show(imgs_ssim, titles=titles_ssim, cols=3, figsize=(14, 14))
The error maps presented were individually normalized (cv2.NORM_MINMAX) to maximize visual contrast and reveal the spatial structure of distortions. This means that:
Therefore, the intensity of white in the maps is NOT comparable across different qualities — the maps serve only to reveal the spatial signature of the error (edges vs. blocks), not its magnitude. The correct magnitude is given by the PSNR and SSIM values, which clearly show that Q=95 has much smaller error than Q=25.
The compression process in the JPEG standard is based on the combined application of spatial, perceptual, and statistical transformations to reduce the redundancies of an image. Table 5.9 summarizes the role of each stage in the pipeline and its respective impact on data reduction.
| Stage | Analytical Operation | Gain / Compression Mechanism |
|---|---|---|
| \(YC_bC_r\) Conversion | Isolation of luminance and chrominance channels. | Models the perception of the HVS, allowing color and brightness to be treated independently. |
| 4:2:0 Subsampling | Reduction of the spatial resolution of the color channels (\(C_b\) and \(C_r\)). | Eliminates approximately 50% of raw data with minimal visual impact. |
| \(8 \times 8\) DCT | Mapping from the spatial domain to the spatial frequency domain. | Energy compaction, concentrating vital information in the first coefficients. |
| Linear Quantization | Integer division of the coefficients by a weighting matrix \(Q(u,v)\). | Main source of lossy compression; eliminates imperceptible high frequencies. |
| Entropy Coding | Application of RLE algorithms and Huffman coding. | Lossless statistical compression, optimized by long runs of null coefficients. |
Applying excessively aggressive compression rates (reduced quality factors) introduces predictable distortions in the reconstructed image, arising from the mathematical limitations of the model:
Bringing together the techniques consolidated throughout this chapter, a complete image restoration pipeline is presented, combining spectral analysis in the frequency domain with adaptive filtering in the spatial domain. The goal is to attenuate mixed noise (composed of Gaussian degradation and periodic interference) while preserving the structural details of the original image as much as possible.
\[ \text{Noisy Image} \xrightarrow{\text{FFT2}} \xrightarrow{\text{Gaussian Notch Filter}} \xrightarrow{\text{IFFT2}} \xrightarrow{\text{Bilateral Filter}} \text{Restored Image} \]
The pair of statistical metrics PSNR and SSIM provides a complementary qualitative and morphological evaluation of the restoration process:
In practice, there is an analytical trade-off between noise reduction and detail preservation: excessively aggressive spatial filters attenuate high-frequency noise well but degrade fine textures and smooth sharp edges — which simultaneously reduces both the PSNR and the SSIM relative to the original image. The challenge of filter design is to find the balance point that maximizes both metrics, ensuring a faithful and visually pleasing restoration.
The numerical and visual results generated by Figure 5.30 demonstrate the practical relevance of combining different processing domains. The simultaneous insertion of periodic and stochastic noise corrupts the morphological properties of the signal, severely reducing the similarity indices and the signal-to-noise ratio of the reference image.
The isolation and suppression of harmonic peaks in the frequency domain through the notch mask remove the sinusoidal interference fringes spread across the two-dimensional space. As evidenced in the printed data of Figure 5.30, this surgical filtering promotes an immediate and substantial leap in the PSNR metric. However, high-frequency Gaussian noise remains uniformly active in the spectrum, requiring a complementary approach.
The final restoration is consolidated in the spatial domain with the introduction of the bilateral filter. Unlike conventional low-pass operators (such as the Gaussian or mean filter), which would indiscriminately smooth noise and structural contours, bilateral filtering computes weights based on geometric proximity and radiometric intensity difference. This adaptive behavior attenuates remaining stochastic fluctuations in smooth transition regions while preserving the sharpness of spatial edges.
The convergence of both approaches results in a substantial and simultaneous improvement in PSNR and SSIM relative to the noisy image — although the final values remain inferior to those of the original image (PSNR = \(\infty\), SSIM = 1.0), due to the inevitable loss of spectral and textural information during the filtering processes. The smooth (Gaussian) attenuation of spectrum peaks avoids ringing artifacts, while the bilateral filter eliminates residual stochastic noise without compromising edge sharpness. The results confirm the effectiveness and practical complementarity of the frequency analysis tools presented in this chapter, demonstrating that hybrid filtering (frequency + spatial) is superior to any isolated approach for restoring images degraded by mixed noise.
import numpy as np
import cv2
try:
from skimage.metrics import structural_similarity as ssim_sk
except ImportError:
import subprocess
subprocess.run(["pip", "install", "scikit-image", "-q"])
from skimage.metrics import structural_similarity as ssim_sk
def suprimir_pico_gaussiano(mask, cy, cx, sigma=3.0):
yy, xx = np.ogrid[:mask.shape[0], :mask.shape[1]]
dist = np.sqrt((yy - cy)**2 + (xx - cx)**2)
notch = np.exp(-dist**2 / (2 * sigma**2))
mask *= (1 - notch)
return mask
# ── 1. Construction of mixed noise ───────────────────────────────────────────
np.random.seed(42)
h_img, w_img = img_gray.shape
X2, Y2 = np.meshgrid(np.arange(w_img), np.arange(h_img))
u0, v0 = 15, 10
ruido_gauss = np.random.normal(0, 15, img_gray.shape)
ruido_period = 30 * np.sin(2 * np.pi * (u0 * X2 / w_img + v0 * Y2 / h_img))
img_noisy = np.clip(img_gray.astype(float) +
ruido_gauss + ruido_period, 0, 255).astype(np.uint8)
# ── 2. Spectrum and identification of peaks ──────────────────────────────────
F_n = np.fft.fftshift(np.fft.fft2(img_noisy.astype(np.float64)))
mag_n = cv2.normalize(np.log1p(np.abs(F_n)), None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)
# ── 3. Gaussian notch on periodic peaks ──────────────────────────────────────
cy0, cx0 = h_img // 2, w_img // 2
mascara_notch = np.ones((h_img, w_img), dtype=np.float64)
for dy, dx in [(+v0, +u0), (-v0, -u0), (+v0, -u0), (-v0, +u0)]:
mascara_notch = suprimir_pico_gaussiano(mascara_notch, cy0 + dy, cx0 + dx, sigma=3.0)
img_notch = np.real(np.fft.ifft2(np.fft.ifftshift(F_n * mascara_notch)))
img_notch = np.clip(img_notch, 0, 255).astype(np.uint8)
# ── 4. Bilateral Filter: removal of residual Gaussian noise ──────────────────
img_den = cv2.bilateralFilter(img_notch, d=7, sigmaColor=25, sigmaSpace=7)
# ── Calculation of Validation Metrics ───────────────────────────────────────
psnr_n, ssim_n = cv2.PSNR(img_gray, img_noisy), ssim_sk(img_gray, img_noisy)
psnr_no, ssim_no = cv2.PSNR(img_gray, img_notch), ssim_sk(img_gray, img_notch)
psnr_d, ssim_d = cv2.PSNR(img_gray, img_den), ssim_sk(img_gray, img_den)
print(f"{'Etapa':>20} | {'PSNR (dB)':>9} | {'SSIM':>6}")
print("-" * 42)
print(f"{'Ruidosa (gauss+per)':>20} | {psnr_n:>9.2f} | {ssim_n:>6.4f}")
print(f"{'Após notch':>20} | {psnr_no:>9.2f} | {ssim_no:>6.4f}")
print(f"{'Notch + bilateral':>20} | {psnr_d:>9.2f} | {ssim_d:>6.4f}")
mascara_vis = (mascara_notch * 255).astype(np.uint8)
mm.show(
[img_gray, img_noisy, mag_n, mascara_vis, img_notch, img_den],
titles=[
"Original",
f"Ruidosa\nPSNR={psnr_n:.1f} dB",
"Espectro\n(picos visíveis)",
"Máscara notch\n(gaussiana suave)",
f"Após notch\nPSNR={psnr_no:.1f} dB",
f"Notch + bilateral\nPSNR={psnr_d:.1f} dB SSIM={ssim_d:.3f}"
],
cols=6, figsize=(20, 4)
) Etapa | PSNR (dB) | SSIM
------------------------------------------
Ruidosa (gauss+per) | 20.19 | 0.3486
Após notch | 24.47 | 0.4715
Notch + bilateral | 28.93 | 0.7426
The transition from the spatial domain to the frequency domain reveals the spectral energy distribution of the image, establishing the analytical foundation for advanced filtering, restoration, and data compression. The structural articulation of these concepts is synthesized in the conceptual map in Figure 5.31.
Next Steps: Chapter 6 inaugurates Part II of the work, applying image processing tools to solve real industrial inspection problems. Techniques for segmentation and shape analysis will be explored for the automatic detection of flaws in production lines—from identifying surface defects on parts to reading QRCode tags on exams, consolidating the bridge between the theory presented in Part I and the practical demands of computer vision.
In this edition, the Gemini Notebook platform is encouraged as a complementary learning tool — not as a substitute for careful reading, problem-solving, or hands-on experimentation. Based on artificial intelligence architectures, the system uses exclusively the instructional material and documents provided by the author as its knowledge base, ensuring that the generated responses are conceptually aligned with the course content and the pedagogical approach adopted throughout this work.
(10%) Direct Implementation of the 2D DFT: Analytically implement the 2D Discrete Fourier Transform (DFT) without the aid of native library functions (such as np.fft.fft2), using strictly the mathematical formulation defined in Equation 5.1 for a matrix of dimensions \(16 \times 16\). Perform numerical validation by comparing the generated coefficients with the results of the np.fft.fft2 function, ensuring that the maximum absolute deviation is less than \(10^{-8}\). Measure the execution times of both methods and provide a theoretical justification for the observed disparity in terms of asymptotic complexity.
(15%) Periodic Noise Suppression: Add sinusoidal interferences with spatial frequencies \((u_0, v_0) \in \{(5,10), (20,5), (30,30)\}\) to the Cameraman test image. For each degradation scenario, design a specific notch filtering mask in the frequency domain to isolate and attenuate the unwanted harmonic peaks. Quantitatively evaluate the effectiveness of the restoration process by computing the PSNR and SSIM metrics. Discuss analytically the trade-off between sinusoidal noise attenuation and the undesirable attenuation of legitimate structural features of the image.
(15%) Comparative Analysis of Low-Pass Operators: Conduct a comparative study among the Ideal, Gaussian, and Butterworth low-pass filters (with harmonic orders \(n = 1, 2, 4\)), parameterized with cutoff frequencies \(D_0 = 20, 40, 60\) pixels. For each structural combination, compute the PSNR and SSIM indices of the resulting image against the original reference signal. Organize the quantitative data in a structured table and plot the one-dimensional graphs of the corresponding transfer functions along the horizontal profile \(H(u, 0)\).
(15%) Haar Multiresolution Filter Bank: Develop a script to manually perform the first-level 2D discrete wavelet decomposition using the Haar family. The algorithm must compute the corresponding low-pass (\(h\)) and high-pass (\(g\)) filter coefficients, applying them separably over the rows and columns of the matrix, followed by the decimation operation (spatial subsampling by a factor of 2). Numerically validate the accuracy of your implementation by comparing the obtained subbands with the output of the pywt.dwt2(img, 'haar') function.
(15%) Sparse Compression by Wavelet Thresholding: Apply the hard thresholding filtering technique to the detail coefficients of the wavelet decomposition, adopting the numerical thresholds \(T \in \{5, 10, 20, 40, 80\}\) for the Haar, Daubechies (db4), and Symlets (sym4) families. After performing the synthesis process through the inverse transform (pywt.waverec2), compute the PSNR and SSIM values of each reconstructed image. Identify and justify which combination of wavelet family and threshold \(T\) maximizes structural similarity.
(15%) Construction of a Simplified JPEG Encoder: Implement the complete data compression pipeline simulating the JPEG standard. The workflow must encompass: spatial conversion \(RGB \rightarrow YC_bC_r\), chroma subsampling at a 4:2:0 ratio, segmentation of luminance into disjoint blocks of \(8 \times 8\) pixels, application of the orthogonal 2D DCT-II, and linear quantization based on the normalized luminance matrix scaled by desired quality factors. Perform the inverse decoding and quantitatively compare the reconstructions with the files generated by the cv2.imencode function for quality factors of 20, 50, and 80.
(15%) Perceptual Analysis on Heterogeneous Content: Develop a synthetic image composed of three distinct regions with contrasting spectral characteristics: a complex photographic texture (representing high stochastic frequencies), a vectorized text area with sharp edges (representing pure step transitions), and a continuous linear gradient (representing homogeneous low frequencies). Submit this mixed image to compression processes under the JPEG, PNG, and WebP formats. Evaluate and interpret the results by correlating the final file size on disk with the obtained PSNR and SSIM metrics, justifying which format exhibits the best performance for signals of a heterogeneous nature and why this advantage occurs in terms of energy compaction and perceptual preservation.
The theoretical foundation and analytical development of the concepts addressed in this chapter are based on the following reference works:
Objective: Apply the concepts presented in Chapter 5 by solving the proposed exercises.
Statement: Write a program that reads an RGB image and converts it to grayscale using the weighted average method. Display the original and converted images side by side.
Hint: Use cv2.cvtColor with cv2.COLOR_BGR2GRAY or implement the formula manually.
Statement: Implement a function to apply a Gaussian blur to an image. Compare the result with the original using different kernel sizes (e.g., 3x3, 5x5, 7x7).
Hint: Use cv2.GaussianBlur with adjustable ksize and sigmaX.
Statement: Perform edge detection using the Canny algorithm. Adjust the thresholds to observe their effect on the detected edges.
Hint: Use cv2.Canny with threshold1 and threshold2 parameters.
Statement: Use morphological operations (erosion and dilation) to remove noise from a binary image. Show the results after each operation.
Hint: Use cv2.erode and cv2.dilate with a structuring element (e.g., cv2.getStructuringElement).
Statement: Develop a simple image segmentation algorithm by thresholding (Otsu’s method). Display the histogram and the resulting binary image.
Hint: Use cv2.threshold with cv2.THRESH_BINARY + cv2.THRESH_OTSU.
After completing all exercises, ensure that the notebook runs without errors. Include comments in your code explaining the main steps. Submit the file according to your course’s guidelines.
End of Chapter 5 Proposed Exercises.
The present list of programming exercises (EP) consolidates the theoretical formulations presented throughout Chapter 5 — Transforms and Compression — through an applied practical track. The exercises are structured around matrices of reduced dimensions, enabling analytical validation and manual inspection of each coefficient, while maintaining the methodological consistency adopted in previous chapters.
The sequencing of the exercises strictly reproduces the conceptual flow of the chapter: it begins with the explicit implementation of the Discrete Fourier Transform (DFT) from its fundamental mathematical definition; it proceeds to the design of low-pass filters and notch masks in the frequency domain; it applies coefficient quantization (the core of lossy compression); and it concludes with the integration of these stages in the construction of a simplified JPEG compression pipeline and in the perceptual analysis of image formats.
In all exercises in this chapter, the coordinates of the spectrum center (the origin of spatial frequencies after the application of the fftshift shift) must be determined using integer division. For a matrix with \(L\) rows and \(C\) columns, the zero-frequency component is located at the position:
\[ (c_y, c_x) = \left( \left\lfloor \frac{L}{2} \right\rfloor, \left\lfloor \frac{C}{2} \right\rfloor \right) \]
This convention is strictly identical to that adopted by the np.fft.fftshift function. Furthermore, in all stages that require discretization or numerical rounding (whether in the quantization of AC coefficients or in the final reconstruction of pixels), the standard rounding to the nearest integer (round half away from zero) must be employed, mitigating ambiguities in values with a fraction exactly equal to \(0.5\).
The notebook enables the development, validation, organization, and testing of solutions for Programming Exercises (PEs) in interactive environments, such as Colab, using the same test cases as Moodle, copying them there only when recording the official grade.
Download morph.py and testsuite.py by running the cell below:
import os, urllib.request
url = "https://raw.githubusercontent.com/fzampirolli/pdi-vc/master/morph/config.py"
if not os.path.exists("config.py"):
urllib.request.urlretrieve(url, "config.py")
import config
config.setup(testsuite=True)
from morph import mm
from testsuite import TestSuite✅ Environment ready. Morph: 1.1.9 | OpenCV: 5.0.0 | TestSuite: 1.1.2To evaluate the tests, run TestSuite("EP05_01.extensão").run() in a new cell, replacing the extension with the one used for your language (.py, .java, .c, .cpp, .js, or .r). The system downloads the test cases from GitHub, runs the program, and calculates the grade automatically.
To test Python code directly, without saving a file, use run_code(codigo) passing the code as a string in a variable named codigo:
codigo = """
from morph import mm
# ... your code here ...
"""
TestSuite("EP05_01").run_code(codigo)In an old document scanner, the sensor captures crumpled paper and fiber texture along with the text — high-frequency noise that “pollutes” the spectrum at the edges. The maintenance technician has no access to the original image, only to the magnitude spectrum already computed by the scanner’s software. Their job is simple and surgical: keep only the central circle of low frequencies (the global structure of the document) and erase everything outside the radius \(D_0\), eliminating the fine texture without even needing to touch the spatial image.
This is the Ideal Low-Pass Filter (LPFI): the most direct spectral operation in the chapter, but also the one that best reveals the anatomy of a centered spectrum.
np.fft.fftshift).| Region | Distance to center | Filter effect |
|---|---|---|
| Center (\(D \le D_0\)) | Low frequencies | Preserved — global structure maintained |
| Edges (\(D > D_0\)) | High frequencies | Zeroed — texture and noise removed |
| Small \(D_0\) | — | Reconstructed image would be very blurry |
| Large \(D_0\) | — | Little filtering; almost all energy preserved |
Input:
Output:
| Input | Output | Observation |
|---|---|---|
| 3 3 1 10 20 30 40 50 60 70 80 90 |
0 20 0 40 50 60 0 80 0 |
Center \((1,1)\). Corners have \(D=\sqrt{2}\approx1.41 > 1\), hence they are zeroed; orthogonal neighbors have \(D=1 \le 1\) and are kept. |
| 1 3 0 5 9 7 |
0 9 0 | \(L=1, C=3\): center at \((0,1)\). Only the central position itself (\(D=0\)) survives \(D_0=0\). |
%%writefile EP05_01.py
# Python codeOverwriting EP05_01.pyTestSuite("EP05_01.py").run()✔️ EP05_01.cases already exists in casos/
📋 5 case(s) loaded from casos/EP05_01.cases
🔍 Testing Python: EP05_01.py
⚠️ EP05_01.py: Empty file (fewer than 3 lines). Tests skipped.
An industrial inspection camera captures images of circuit boards, but the production line’s power supply introduces a periodic electrical interference — a stripe pattern almost imperceptible to the naked eye, yet visible in the Fourier spectrum as pairs of bright peaks symmetrically positioned around the center. The computer vision team cannot redo the capture: they must surgically locate and erase these peak pairs in the spectrum, preserving all other useful image information.
This is the role of the notch reject filter: unlike a low-pass filter (which affects a continuous region), it targets specific points and their symmetric counterparts, leaving the rest of the spectrum untouched.
| Concept | Role in the notch filter |
|---|---|
| Peak at \((\Delta v, \Delta u)\) | Frequency of the periodic interference visually detected in the spectrum |
| Symmetric point \((-\Delta v,-\Delta u)\) | Every DFT of a real signal is Hermitian: peaks always appear in pairs symmetric about the center |
| Radius \(r\) | Controls the “width” of rejection — a large \(r\) removes more energy around the peak, but also useful information |
Input:
Output:
| Input | Output | Observation |
|---|---|---|
| 5 5 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 1 1 1 0 |
1 2 3 4 5 6 0 8 9 10 11 12 13 14 15 16 17 18 0 20 21 22 23 24 25 |
Center \((c_y, c_x) = (2, 2)\). Reported peak \((\Delta v, \Delta u) = (1, 1)\) generates the point \((3, 3)\) (value 19) and its symmetric counterpart \((1, 1)\) (value 7), both zeroed with \(r=0\) (only the exact points). |
%%writefile EP05_02.py
# Python codeOverwriting EP05_02.pyTestSuite("EP05_02.py").run()✔️ EP05_02.cases already exists in casos/
📋 5 case(s) loaded from casos/EP05_02.cases
🔍 Testing Python: EP05_02.py
⚠️ EP05_02.py: Empty file (fewer than 3 lines). Tests skipped.
A photo gallery application needs to reduce the size of thousands of images before uploading them to the cloud, without recoding everything from scratch. The engineer in charge already has the DCT coefficients for each \(4\times4\) block calculated (the computationally expensive step has already been done) — only the quantization table needs to be applied, the step that actually discards information and generates compression. High-frequency coefficients, which are less perceptible to the human eye, receive large divisors and tend to become zero; low-frequency coefficients, which are more perceptible, receive small divisors and survive almost intact.
You will implement exactly this step: quantize and dequantize (divide, round, multiply back) — the heart of JPEG lossy compression.
.5 values never occur in the test cases).| Coefficient | Frequency | Typical value of \(Q\) | Effect of quantization |
|---|---|---|---|
| \(C(0,0)\) | DC (block average) | Small | Almost always survives — dominates energy |
| \(C(u,v)\) low \(u+v\) | Low frequency | Small/medium | Partially preserved |
| \(C(u,v)\) high \(u+v\) | High frequency | Large | Often becomes zero — source of compression |
Input:
Output:
| Input | Output | Observation |
|---|---|---|
| 4 50 10 -5 0 8 -3 2 1 0 1 0 0 2 0 0 -1 2 5 7 8 4 7 8 11 6 8 11 12 9 11 12 14 |
50 10 -7 0 8 0 0 0 0 0 0 0 0 0 0 0 |
\(C(0,0)=50/2=25 \to 25\times2=50\) (preserved). \(C(0,2)=-5/7\approx-0.71\to-1\to-1\times7=-7\). Meanwhile, \(C(1,1)=-3/7\approx-0.43\to0\): zeroed by quantization — most of the block becomes zero, illustrating the energy compaction in the upper-left corner. |
%%writefile EP05_03.py
# Python codeOverwriting EP05_03.pyTestSuite("EP05_03.py").run()✔️ EP05_03.cases already exists in casos/
📋 5 case(s) loaded from casos/EP05_03.cases
🔍 Testing Python: EP05_03.py
⚠️ EP05_03.py: Empty file (fewer than 3 lines). Tests skipped.
A research laboratory in computational astronomy received, from an old mission, a small experimental sensor whose raw data cannot be processed by modern FFT libraries — the validation environment is isolated and only allows basic arithmetic operations. The team needs to reimplement the 2D Discrete Fourier Transform from the mathematical definition itself, cell by cell, to later compare bit by bit with np.fft.fft2 in another environment.
This is the most conceptual exercise on the list: there are no shortcuts. You will implement the double summation from Equation 5.1 directly, demonstrating why the FFT exists — and the computational cost it avoids.
fftshift — the DC remains at \((0,0)\)).fftshift: the output maintains the raw DFT convention, with the DC component at \(F(0,0)\) (upper-left corner)..5 ambiguity.| Element | Meaning |
|---|---|
| \(F(0,0)\) | DC component — sum of all pixels, \(F(0,0) = \sum f(x,y)\) |
| Real part \(\text{Re}(F)\) | Projection of the signal onto cosines |
| Imaginary part \(\text{Im}(F)\) | Projection of the signal onto sines |
| Complexity of this implementation | \(\mathcal{O}((MN)^2)\) — this is why the FFT, with \(\mathcal{O}(MN\log(MN))\), is indispensable for real images |
Input:
Output:
| Input | Output | Observation |
|---|---|---|
| 2 2 1 2 3 4 |
10 2 4 0 |
\(F(0,0)=1+2+3+4=10\) (DC = total sum). \(F(0,1)=(1-2)+(3-4)=-2 \to |F|=2\). \(F(1,0)=(1+2)-(3+4)=-4\to|F|=4\). \(F(1,1)=(1-2)-(3-4)=0\). |
%%writefile EP05_04.py
# Python codeOverwriting EP05_04.pyTestSuite("EP05_04.py").run()✔️ EP05_04.cases already exists in casos/
📋 5 case(s) loaded from casos/EP05_04.cases
🔍 Testing Python: EP05_04.py
⚠️ EP05_04.py: Empty file (fewer than 3 lines). Tests skipped.
You have been hired to create, from scratch, a didactic JPEG codec in an embedded environment, without any available image library—only basic mathematical operations. The client wants to understand exactly where quality is lost and where it is recovered, block by block. This is the final challenge of the chapter: integrate everything that has been studied—the orthonormal DCT-II, perceptual quantization, and IDCT-based reconstruction—into a single end-to-end pipeline, processing an \(N \times N\) block from start to finish, exactly as the JPEG standard does internally, \(8\times8\) pixels at a time.
.5 ambiguity.| Stage | Analogous in the real JPEG standard | Where quality is lost |
|---|---|---|
| Centering | Same—DCT assumes a signal centered at zero | No loss |
| DCT-II | Steps 3–4 of the pipeline (Table 5.7) | No loss (exact and reversible transformation) |
| Quantization | Step 5—division by \(Q(u,v)\) | Main source of loss—high-frequency coefficients become zero |
| IDCT | Final reconstruction | Reconstructs exactly the quantized coefficients, not the original ones |
Input:
Output:
| Input | Output | Observation |
|---|---|---|
| 4 120 130 125 128 115 140 135 122 118 150 160 130 110 120 145 138 4 6 8 10 6 8 10 12 8 10 12 16 10 12 16 20 |
118 126 119 131 114 143 140 119 117 149 159 130 107 121 146 139 |
After DCT, aggressive quantization of high frequencies (large \(Q\) values in the lower-right corner) and IDCT-based reconstruction, the block remains close to the original but not identical—the difference is the cost of lossy compression. |
If the result does not match, check in this order: (1) the raw DCT coefficients (before quantization)—they should reconstruct the original exactly via IDCT if you skip steps 6–7; (2) the \(\alpha(u)\) table—a common mistake is applying \(\sqrt{2/N}\) also for \(u=0\); (3) the quantization rounding, which must occur before multiplying back by \(Q\).
%%writefile EP05_05.py
# Python codeOverwriting EP05_05.pyTestSuite("EP05_05.py").run()✔️ EP05_05.cases already exists in casos/
📋 5 case(s) loaded from casos/EP05_05.cases
🔍 Testing Python: EP05_05.py
⚠️ EP05_05.py: Empty file (fewer than 3 lines). Tests skipped.