Nei capitoli precedenti, tutte le operazioni sono state eseguite nel dominio spaziale, in cui gli algoritmi agiscono direttamente sui valori di intensità dei pixel.
In questo capitolo verrà presentato un approccio complementare: il dominio della frequenza, nel quale l’immagine è rappresentata dalle variazioni spaziali di intensità, e non solo dai valori individuali dei pixel.
Il concetto di frequenza spaziale descrive la rapidità con cui l’intensità varia lungo l’immagine. Le variazioni lente corrispondono a basse frequenze, mentre bordi, dettagli fini e rumori corrispondono ad alte frequenze.
Questa rappresentazione si basa sul fatto che qualsiasi immagine digitale discreta può essere decomposta in una combinazione di funzioni ortogonali. La Trasformata di Fourier utilizza una base di esponenziali complesse bidimensionali (equivalenti a sinusoidi con orientamento e frequenza specifici). Altre trasformate, come la Trasformata del Coseno (DCT) e la Trasformata Wavelet (DWT), utilizzano diverse famiglie di funzioni di base — coseni bidimensionali nel caso della DCT, e funzioni con supporto compatto nel caso delle wavelet.
Tra le principali applicazioni di questa rappresentazione si evidenziano:
Al termine di questo capitolo, sarai in grado di:
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✅ Ambiente pronto. Morph: 1.1.9 | OpenCV: 5.0.0L’analisi di Fourier si basa sul principio secondo cui qualsiasi segnale periodico può essere rappresentato come una somma di funzioni sinusoidali con differenti frequenze, ampiezze e fasi. Questo concetto si applica anche alle immagini digitali, consentendo di rappresentarle nel dominio della frequenza invece che nel dominio spaziale.
La Figura 5.1 illustra questa decomposizione per un segnale unidimensionale. Nel caso di un’immagine, la Trasformata Discreta di Fourier (DFT) converte la matrice delle intensità \(f(x,y)\) in un insieme di coefficienti che descrive il contributo delle diverse frequenze spaziali presenti nell’immagine.
Prima di studiare le immagini bidimensionali, il simulatore della Figura 5.2 illustra il principio dell’analisi di Fourier per segnali unidimensionali: una forma d’onda può essere approssimata dalla somma di sinusoidi con diverse frequenze e ampiezze.
Man mano che si aggiungono nuovi termini, la somma delle sinusoidi (curva nera) si avvicina alla forma d’onda di riferimento (tratteggiata). Il grafico inferiore mostra lo spettro delle ampiezze, indicando il contributo di ciascuna frequenza alla ricostruzione del segnale.
Esplora il simulatore e rispondi:
1. Quanti termini sono necessari per una buona approssimazione dell’onda quadra?
Con circa 15-20 termini, la forma dell’onda si avvicina già bene al riferimento. Tuttavia, in prossimità delle discontinuità permane una piccola oscillazione, nota come fenomeno di Gibbs, che non scompare nemmeno aggiungendo più termini.
2. Quale forma converge più rapidamente? Perché?
L’onda triangolare converge più rapidamente, poiché le ampiezze delle sue armoniche decadono più velocemente di quelle dell’onda quadra e dell’onda a dente di sega. Di conseguenza, pochi termini già producono una buona approssimazione.
3. Come cambia lo spettro tra l’onda quadra e quella triangolare?
Entrambe presentano solo armoniche dispari, ma nell’onda triangolare le ampiezze diminuiscono molto più rapidamente. Pertanto, poche armoniche sono sufficienti per ricostruire il segnale con buona precisione.
Applicando la Trasformata Discreta di Fourier (DFT) a un’immagine e visualizzando il modulo dei suoi coefficienti (vedi Figura 5.5), si ottiene lo spettro di ampiezza, che mostra la distribuzione delle frequenze spaziali presenti nell’immagine.
Il coefficiente situato all’origine della DFT, denominato componente DC (Direct Current), corrisponde alla frequenza nulla e rappresenta l’intensità media dell’immagine. Per convenzione, tale coefficiente è memorizzato nell’angolo superiore sinistro dello spettro. Per facilitarne l’interpretazione, si applica l’operazione FFT Shift, che sposta la componente DC al centro dell’immagine. Dopo questo spostamento, le basse frequenze si concentrano nella regione centrale, mentre le alte frequenze si trovano in prossimità dei bordi, come riassunto nella Tabella 5.1.
| Regione dello spettro | Componenti predominanti | Esempi nell’immagine |
|---|---|---|
| Centro (basse frequenze) | Variazioni spaziali lente | Illuminazione, regioni omogenee e forme globali |
| Regione intermedia (medie frequenze) | Variazioni di scala intermedia | Trame e pattern ripetitivi |
| Bordi (alte frequenze) | Variazioni spaziali rapide | Contorni, dettagli fini e rumore |
Questa organizzazione facilita l’interpretazione dello spettro e la progettazione di filtri. L’attenuazione delle basse frequenze riduce le variazioni globali di intensità, mentre l’attenuazione delle alte frequenze smussa l’immagine riducendo i dettagli fini e parte del rumore.
Prima di presentare la formulazione matematica della Trasformata Discreta di Fourier (DFT), è utile analizzare la sua inversa, denominata Trasformata Discreta Inversa di Fourier (IDFT). Si consideri uno spettro in cui tutti i coefficienti siano nulli, eccetto uno. Un esempio di questa costruzione è presentato nel codice della Figura 5.3 e può essere esplorato interattivamente nel simulatore della Figura 5.4..
L’immagine ricostruita è una sinusoide bidimensionale. La posizione del coefficiente nello spettro determina la sua orientazione e la sua frequenza spaziale, mentre la sua magnitudine e la sua fase definiscono, rispettivamente, la sua ampiezza e il suo spostamento spaziale. Pertanto, ciascun coefficiente della DFT rappresenta una componente sinusoidale, e l’immagine originale può essere ricostruita mediante la somma di tutte queste componenti.
import cv2
N_grid = 100
espectro_vazio = np.zeros((N_grid, N_grid), dtype=complex)
# Accendendo un singolo punto (frequenza) fuori dal centro
u0, v0 = 10, 5
espectro_vazio[N_grid//2 - v0, N_grid//2 - u0] = 1000
# Tornando al dominio spaziale (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)
# Evidenziazione visiva del punto
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=["Spettro (1 punto attivo)", "Onda 2D risultante (IDFT)"],
cols=2, figsize=(10, 4))
Si consideri un’immagine \(f(x,y)\) con dimensioni \(M \times N\). La sua Trasformata Discreta di Fourier 2D (DFT) è definita da:
\[ 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}\]
dove \(u = 0, 1, \ldots, M-1\) e \(v = 0, 1, \ldots, N-1\) rappresentano le frequenze discrete nelle direzioni orizzontale e verticale, rispettivamente. Il termine esponenziale corrisponde a una sinusoide bidimensionale, la cui frequenza e orientazione sono determinate dagli indici \((u,v)\).
La Trasformata Discreta Inversa di Fourier 2D (IDFT) ricostruisce l’immagine originale a partire dai suoi coefficienti:
\[ 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}\]
Le Equazioni Equazione 5.1 e Equazione 5.2 mostrano che la DFT e la IDFT formano una coppia di trasformazioni: la prima converte l’immagine nel dominio della frequenza, mentre la seconda ricostruisce esattamente l’immagine originale a partire dai suoi coefficienti.
Il termine \(j\) denota l’unità immaginaria, definita da \(j^2 = -1\). In ingegneria e nell’elaborazione dei segnali, si adotta \(j\) invece di \(i\) per evitare conflitti con la notazione della corrente elettrica. Il suo utilizzo nell’esponenziale complessa, governato dalla formula di Eulero (\(e^{j\theta} = \cos\theta + j\sin\theta\)), consente di rappresentare in modo compatto l’ampiezza e la fase di ciascuna frequenza spaziale presente nell’immagine.
Il coefficiente \(F(0,0)\), denominato componente DC (Direct Current), è uguale alla somma delle intensità di tutti i pixel dell’immagine (vedi Figura 5.5):
\[ F(0,0)=MN\,\bar{f}, \]
dove \(\bar{f}\) è l’intensità media dell’immagine. Per questo motivo, il componente DC rappresenta il livello medio di intensità e, nella maggior parte delle immagini naturali, possiede la maggiore ampiezza dello spettro.
Gli altri coefficienti rappresentano variazioni attorno a tale media. Dopo l’applicazione del FFT Shift, il componente DC viene spostato al centro dello spettro, concentrando le basse frequenze nella regione centrale e le alte frequenze ai bordi.
Ogni coefficiente della Trasformata Discreta di Fourier (DFT) è un numero complesso e può essere scritto come
\[ F(u,v)=R(u,v)+j\,I(u,v), \]
dove \(R(u,v)\) e \(I(u,v)\) corrispondono, rispettivamente, alle parti reale e immaginaria del coefficiente. Dall’Equazione Equazione 5.1 si ottiene
\[ 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), \]
e
\[ 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). \]
Da questa rappresentazione si definiscono due grandezze fondamentali:
\[ |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). \]
Pertanto, ogni coefficiente può anche essere scritto nella sua forma polare,
\[ F(u,v)=|F(u,v)|\,e^{j\phi(u,v)}. \]
Lo spettro di Fourier può quindi essere visualizzato mediante due immagini distinte: lo spettro di magnitudine, generalmente utilizzato per analizzare la distribuzione delle frequenze, e lo spettro di fase, che descrive l’organizzazione spaziale delle componenti sinusoidali.
Sebbene lo spettro di magnitudine sia il più utilizzato per l’ispezione visiva, la fase contiene gran parte delle informazioni strutturali dell’immagine. La combinazione di magnitudine e fase consente di ricostruire esattamente l’immagine originale tramite la IDFT.
Una dimostrazione classica consiste nel combinare l’ampiezza di un’immagine con la fase di un’altra e ricostruire il risultato. Questo esperimento evidenzia che:
Quando un’immagine viene ricostruita con l’ampiezza di A e la fase di B, il risultato tende a somigliare più a B che ad A, evidenziando che la fase è il componente principale responsabile dell’organizzazione spaziale della scena. Tuttavia, l’ampiezza rimane importante, poiché modula il contrasto delle strutture ricostruite. Pertanto, una ricostruzione fedele dipende dalla combinazione coerente tra ampiezza e fase.
Un esempio di questo comportamento è presentato in Figura 5.6.
La fase di un segnale svolge ruoli distinti nell’audio e nelle immagini:
In entrambi i domini, l’ampiezza è legata all’intensità delle componenti di frequenza: nell’audio, influenza il timbro e l’intensità percepita; nelle immagini, influenza il contrasto e la tessitura.
# ── Esperimento: L'Importanza della Fase ─────────────────────────────────────
# ── Caricamento dell'immagine ────────────────────────────────────────────────
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))
# Creare un'immagine B sintetica (motivo geometrico)
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)
# Scambio di Fase
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=["Immagine A", "Immagine B", "Mag(A) + Fase(B)", "Mag(B) + Fase(A)"],
cols=4, figsize=(16, 4)
)
print("💡 La fase preserva bordi e contorni; la magnitudine controlla contrasto e")
print("texture. Nell'audio stereo, la fase influisce sulla localizzazione spaziale; nelle")
print("immagini, determina l'organizzazione della scena.")
💡 La fase preserva bordi e contorni; la magnitudine controlla contrasto e
texture. Nell'audio stereo, la fase influisce sulla localizzazione spaziale; nelle
immagini, determina l'organizzazione della scena.Il Teorema della Convoluzione stabilisce una relazione fondamentale tra il dominio spaziale e quello delle frequenze:
\[ f(x,y) \circledast h(x,y) \;\overset{\mathcal{F}}{\longleftrightarrow}\; F(u,v)\,H(u,v) \tag{5.3}\]
dove \(\circledast\) rappresenta la convoluzione circolare discreta. Pertanto, la convoluzione tra un’immagine \(f(x,y)\) e un filtro \(h(x,y)\) può essere sostituita dalla moltiplicazione dei loro spettri.
In pratica, per ottenere lo stesso risultato della convoluzione lineare eseguita nel dominio spaziale, si applica il zero-padding prima della Trasformata Rapida di Fourier (FFT), evitando artefatti ai bordi dell’immagine.
Tuttavia, la filtrazione nel dominio delle frequenze non è sempre l’alternativa più efficiente. Per filtri come quello gaussiano e il filtro della media (Box Filter), la proprietà di separabilità consente di ridurre significativamente il costo computazionale della convoluzione nel dominio spaziale.
Un kernel separabile può essere scritto come il prodotto esterno di due vettori unidimensionali,
\[ H = v\,h^T, \]
consentendo di sostituire la convoluzione bidimensionale con due convoluzioni unidimensionali consecutive: una nella direzione orizzontale e una in quella verticale.
Un kernel non separabile, invece, non ammette tale decomposizione e, pertanto, la sua convoluzione deve essere eseguita direttamente sull’intorno bidimensionale.
In pratica, per un kernel di dimensione \(K \times K\), la convoluzione diretta richiede \(K^2\) moltiplicazioni per pixel, mentre un kernel separabile ne richiede solo \(2K\), riducendo significativamente il costo computazionale.
Si consideri un’immagine di dimensioni \(M \times N\) e un filtro quadrato di dimensione \(K \times K\). La Tabella 5.2 confronta la complessità delle principali strategie di filtraggio.
| Metodo di filtraggio | Complessità asintotica | Dipendenza da \(K\) | Applicazione tipica |
|---|---|---|---|
| Spatiale non separabile | \(\mathcal{O}(MNK^2)\) | Quadratica | Kernel piccoli e non separabili |
| Spatiale separabile | \(\mathcal{O}(MNK)\) | Lineare | Filtri Gaussiano e della media |
| Tramite FFT | \(\mathcal{O}(MN\log(MN))\) | Indipendente da \(K\) | Kernel grandi |
Per kernel piccoli, la convoluzione spaziale, specialmente quando il filtro è separabile, risulta solitamente più efficiente grazie al basso costo delle operazioni. All’aumentare della dimensione del kernel, il filtraggio tramite FFT diventa più vantaggioso, poiché il suo costo è praticamente indipendente dalla dimensione del filtro.
Il grafico ottenuto nella prova con l’immagine delle monete (\(2560 \times 1920\)), presentato nella Figura 5.7, conferma il comportamento previsto dall’analisi della complessità computazionale.
Convoluzione non separabile (\(\mathcal{O}(MNK^2)\))
La convoluzione diretta presenta una crescita quadratica con la dimensione del kernel. Per valori piccoli di \(K\), il costo è basso, ma aumenta rapidamente man mano che il kernel cresce, diventando impraticabile per applicazioni in tempo reale.
Filtraggio tramite FFT (\(\mathcal{O}(MN \log(MN))\))
Il costo della FFT dipende solo dalla dimensione dell’immagine, essendo indipendente da \(K\). Pertanto, le sue prestazioni rimangono approssimativamente costanti al variare del kernel, rendendola vantaggiosa per filtri grandi o non separabili.
Convoluzione separabile (\(\mathcal{O}(MNK)\))
La scomposizione del kernel in due filtri monodimensionali riduce significativamente il costo computazionale. In pratica, questo approccio tende ad essere il più efficiente per filtri separabili, specialmente in implementazioni ottimizzate.
In generale, la scelta del metodo dipende dalla dimensione e dalla struttura del kernel. I filtri separabili sono più efficienti nel dominio spaziale, mentre la FFT diventa più vantaggiosa per kernel grandi o per convoluzioni multiple nel dominio della frequenza.
\[ g = \mathcal{F}^{-1}\bigl[\mathcal{F}(f)\cdot \mathcal{F}(h)\bigr] \quad \text{(FFT)} \qquad g = f \circledast h \quad \text{(convoluzione diretta)} \qquad g = (f \circledast v) \circledast h^T \quad \text{(separabile)} \tag{5.4}\]
dove:
La Trasformata Discreta di Fourier (DFT) assume che l’immagine sia estesa periodicamente nello spazio, cioè che i suoi bordi si ripetano indefinitamente.
In questa condizione, la moltiplicazione nel dominio della frequenza corrisponde a una convoluzione circolare nel dominio spaziale. Di conseguenza, regioni opposte dell’immagine (alto e basso, sinistra e destra) iniziano a interagire artificialmente, come illustrato nella Figura 5.8..
L’applicazione del zero-padding prima della FFT riduce questo effetto estendendo l’immagine con valori nulli ai bordi, avvicinando il risultato alla convoluzione lineare. Questo comportamento può essere interpretato alla luce del Teorema della Convoluzione, presentato nella Figura 5.9..
# ── Definire M e N ─────────────────────────────────────────────────────────────
M, N = img_gray.shape # riga aggiunta
# Simulazione di un filtro di dislocamento brutale
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))
# Filtraggio SENZA padding (causa il 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=["Originale", "Filtraggio s/ Padding (Perdita)"], cols=2, figsize=(10, 4))
# ── Kernel Gaussiano 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) # kernel 2D separabile
# ── Metodo 1: Convoluzione spaziale diretta ─────────────────────────────────────
f_float = img_gray.astype(np.float64)
conv_esp = cv2.filter2D(f_float, -1, kernel, borderType=cv2.BORDER_CONSTANT)
# ── Metodo 2: Moltiplicazione in frequenza (via FFT) ──────────────────────────
M, N = f_float.shape
# Padding per convoluzione lineare (evita aliasing circolare)
Mpad = 2 ** int(np.ceil(np.log2(M + K - 1)))
Npad = 2 ** int(np.ceil(np.log2(N + K - 1)))
# Posiziona il kernel con l'origine in (0,0) e padding con zeri
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))
# Ritaglio per compensare lo spostamento introdotto dal posizionamento del kernel
offset = K // 2
conv_freq_crop = conv_freq[offset:offset+M, offset:offset+N]
# ── Verifica numerica ──────────────────────────────────────────────────────
diff = np.abs(conv_esp - conv_freq_crop)
print(f"Differenza massima (|conv_spaz - conv_freq|): {diff.max():.2e}")
print(f"Differenza media (|conv_spaz - conv_freq|): {diff.mean():.2e}")
print(f"→ Teorema della Convoluzione verificato numericamente.")
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=[
"Originale",
"Convoluzione spaziale",
"Moltiplicazione in frequenza",
f"Differenza (máx={diff.max():.1e})"
],
cols=4, figsize=(16, 5)
)Differenza massima (|conv_spaz - conv_freq|): 2.56e-13
Differenza media (|conv_spaz - conv_freq|): 2.76e-14
→ Teorema della Convoluzione verificato numericamente.
La differenza residua dell’ordine di \(10^{-13}\) non viola il Teorema della Convoluzione, ma riflette limitazioni computazionali inerenti all’aritmetica a virgola mobile (doppia precisione, ~\(10^{-16}\)) e all’ordine delle operazioni tra i due metodi:
Pertanto, l’uguaglianza teorica è esatta, ma l’implementazione numerica produce una differenza praticamente nulla (errore relativo < \(10^{-12}\)), confermando il teorema entro la precisione della macchina.
Un filtro nel dominio delle frequenze può essere interpretato come una funzione di trasferimento applicata allo spettro dell’immagine. In questa rappresentazione, ogni coefficiente di frequenza viene moltiplicato per un valore compreso tra 0 e 1, che determina la sua attenuazione o preservazione. La forma di questa funzione definisce l’effetto visivo del filtro.
Taglio netto e ringing. I filtri ideali con transizione istantanea a una frequenza di taglio \(D_0\) producono discontinuità nel dominio delle frequenze. Questa discontinuità si riflette nel dominio spaziale come oscillazioni in prossimità dei bordi, note come ringing. Questo effetto è associato alla convoluzione con funzioni di supporto infinito nello spazio, come la funzione sinc, come illustrato nella Figura 5.10.
Filtri con transizione graduale. Alternative come i filtri Gaussiano e Butterworth attenuano la transizione tra regioni di passaggio e di reiezione, riducendo il ringing. Di contro, tale attenuazione implica un confine di separazione meno definito tra frequenze preservate e attenuate.
# Simulando il Filtro Ideale in Frequenza (Cilindro) e la sua rappresentazione Spaziale (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)
# Frequenza: Cilindro Ideale (1 al centro, 0 fuori dal raggio 20)
H_freq = np.zeros((N_grid, N_grid))
H_freq[D <= 20] = 1
# Spazio: L'inversa risulta nella famigerata Sinc 2D
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()
I filtri passa-basso attenuano le componenti ad alta frequenza, producendo un appiattimento dell’immagine e una riduzione del rumore. Dopo la centralizzazione dello spettro (FFT Shift), la distanza di ciascun punto dal centro è data da:
\[ D(u,v) = \sqrt{\left(u - \tfrac{M}{2}\right)^2 + \left(v - \tfrac{N}{2}\right)^2} \tag{5.5}\]
Filtro Ideale (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}\]
Il taglio netto a \(D_0\) introduce discontinuità nel dominio delle frequenze, causando oscillazioni nel dominio spaziale note come ringing. Questo effetto è associato alla convoluzione con funzioni a supporto infinito.
Filtro Gaussiano (LPFG): \[ H_{\text{gauss}}(u,v) = e^{-D^2(u,v)/(2\sigma^2)} \tag{5.7}\]
La regolarità della funzione gaussiana nel dominio delle frequenze evita discontinuità, eliminando il ringing e producendo una transizione graduale tra frequenze preservate e attenuate.
Filtro di Butterworth (LPFB) di ordine \(n\): \[ H_{\text{BW}}(u,v) = \frac{1}{1 + \left[D(u,v)/D_0\right]^{2n}} \tag{5.8}\]
Il parametro \(n\) controlla la gradualità della transizione tra l’attenuazione e il passaggio delle frequenze. Valori piccoli producono transizioni morbide, mentre valori grandi avvicinano il comportamento al filtro ideale, con un maggiore rischio di ringing. Un esempio comparativo è mostrato in Figura 5.11.
I filtri passa-alto possono essere ottenuti da un filtro passa-basso complementare, definito come:
\[ H_{\text{HP}}(u,v) = 1 - H_{\text{LP}}(u,v) \]
Questo tipo di filtro preserva le componenti ad alta frequenza, evidenziando bordi e dettagli, mentre attenua le regioni a variazione graduale.
I filtri passa-banda preservano solo una banda intermedia di frequenze, limitata da due raggi \(D_L\) e \(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] \]
Questo tipo di filtraggio è utile quando si desidera rimuovere simultaneamente le componenti a bassa e ad alta frequenza, preservando solo le strutture di scala intermedia.
Un’applicazione importante è la rimozione del rumore periodico, in cui i pattern regolari appaiono come picchi localizzati nello spettro di magnitudine. Questi picchi possono essere attenuati mediante filtri notch (reietta-banda), posizionati specificamente sulle frequenze indesiderate.
Esempi di filtri nel dominio della frequenza sono presentati nel simulatore della Figura 5.12, Figura 5.13 e Figura 5.14..
import io
def distancia_centro(M, N):
"""Matrice delle distanze dal centro dello spettro."""
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):
"""Applica il filtro H (centrato) all'immagine tramite 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 # frequenza di taglio
n_bw = 2 # ordine del Butterworth
# ── Funzioni di trasferimento ──────────────────────────────────────────────────
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))
# ── Immagini filtrate ─────────────────────────────────────────────────────────
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)
# ── Profili di H(u,v) ─────────────────────────────────────────────────────────
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)
# ── Filtri H visualizzati ────────────────────────────────────────────────────
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=[
"Originale", "LPF Ideale", "LPF Gaussiano", "LPF Butterworth (n=2)",
"H Ideale", "H Gaussiano","H Butterworth", "Profili H(u,v)"
],
cols=4, figsize=(16, 9)
)
# Filtro passa-alta: complemento do passa-baixa Gaussiano
# Reutiliza aplicar_filtro_freq() definida na célula anterior
H_alta = 1 - H_gauss
img_alta = aplicar_filtro_freq(img_gray, H_alta)
mm.show(
[img_gray, img_alta],
titles=["Originale", "Passa-alto Gaussiano ($D_0=30$)"],
cols=2
)
Il rumore periodico — associato a interferenze elettriche, pattern regolari dei sensori o artefatti di scansione — appare nello spettro di Fourier come picchi puntuali simmetrici rispetto al centro.
Il filtro reietta-banda (notch) attenua selettivamente queste frequenze, preservando le altre componenti dell’immagine. Un esempio di applicazione è presentato nella Figura 5.15..
# ── Immagine con rumore periodico sintetico ──────────────────────────────────────
h_img, w_img = img_gray.shape
x = np.arange(w_img)
y = np.arange(h_img)
X, Y = np.meshgrid(x, y)
# Usare frequenze intere e consistenti (importante per una restaurazione perfetta)
u0, v0 = 20, 20 # frequenze esatte del rumore
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)
# ── Spettro dell'immagine rumorosa ───────────────────────────────────────────────
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)
# ── Maschera notch ────────────────────────────────────────────────────────────
mascara = np.ones((h_img, w_img), dtype=np.float64)
r_notch = 8 # raggio del notch (regolazione fine se necessario)
def suprimir_pico(mask, cy, cx, r):
"""Azera un disco di raggio r centrato in (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
# Coordinate centrali
cy, cx = h_img // 2, w_img // 2
# Sopprimere i 4 picchi simmetrici (importante!)
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)
# ── Filtraggio e ricostruzione ─────────────────────────────────────────────────
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)
# Valutazione
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 (originale vs restaurata): {psnr:.2f} dB")
# ── Visualizzazione ─────────────────────────────────────────────────────────────
mm.show(
[img_ruidosa, mag_vis, mascara_vis, img_rest_vis],
titles=[
"Con rumore periodico",
"Spettro (log)",
"Maschera notch",
f"Restaurata (PSNR={psnr:.1f} dB)"
],
cols=4,
figsize=(16, 4)
)PSNR (originale vs restaurata): 32.93 dB
| Filtro | Effetto visivo | Artefatto | Uso |
|---|---|---|---|
| Passa-basso ideale | Attenuazione intensa | Ringing | Illustrativo |
| Passa-basso Gaussiano | Attenuazione graduale | Non presenta ringing | Attenuazione generale |
| Passa-basso Butterworth | Attenuazione controllata | Ringing (ordini elevati) | Compromesso tra attenuazione e selettività |
| Passa-alto | Enfatizzazione dei bordi | Amplificazione del rumore | Rilevamento dei contorni |
| Notch | Rimozione selettiva delle frequenze | Possibili distorsioni locali | Rimozione del rumore periodico |
La progettazione dei filtri nel dominio della frequenza consiste nella definizione di maschere spettrali. Tuttavia, effetti nel dominio spaziale, come ringing e sfocatura, emergono direttamente da queste scelte nello spettro.
La Trasformata di Fourier decompone il segnale in frequenze globali: ogni coefficiente \(F(u,v)\) riceve contributi dall’intera immagine, senza informazioni esplicite sulla localizzazione spaziale di tali frequenze. Pertanto, strutture localizzate, come i bordi, sono rappresentate in modo distribuito nello spettro.
Le wavelet (ondine) superano questa limitazione utilizzando funzioni base localizzate nello spazio, che possono essere traslate e scalate. Queste funzioni possiedono supporto compatto, ossia sono diverse da zero solo in una regione finita del dominio, consentendo una rappresentazione simultanea in termini di frequenza e localizzazione spaziale.
La Trasformata di Fourier descrive con precisione quali frequenze sono presenti in un segnale, ma non rappresenta esplicitamente dove queste frequenze si verificano nello spazio.
Nell’esperimento presentato nella Figura 5.16, due immagini con strutture localizzate in posizioni diverse producono spettri di magnitudine praticamente identici. Ciò accade perché la rappresentazione di Fourier è globale: ogni coefficiente riceve un contributo dall’intera immagine.
Di conseguenza, lo spettro di magnitudine non rappresenta esplicitamente la localizzazione dei bordi o di altre strutture, ma solo la distribuzione delle frequenze presenti. Questa limitazione ha motivato lo sviluppo di rappresentazioni multirisoluzione, come la Trasformata Wavelet Discreta (DWT), in grado di descrivere simultaneamente la frequenza e la localizzazione spaziale delle strutture dell’immagine.
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=["Señale A", "Spettro A", "Señale B (Spostato)", "Spettro B"],
cols=4, figsize=(14, 4)
)
La Trasformata Wavelet Discreta (DWT) applica, separatamente nelle direzioni orizzontale e verticale, due filtri complementari: un passa-basso \(h\) (approssimazione) e un passa-alto \(g\) (dettagli), seguiti da sottocampionamento per un fattore di 2 in ciascuna dimensione. Questo processo produce quattro sottobande, i cui nomi indicano la combinazione dei filtri applicati in ciascuna direzione (L = Low-pass, passa-basso; H = High-pass, passa-alto). Le caratteristiche di ciascuna sottobanda sono riassunte nella Tabella 5.3.
\[ \text{DWT}(f)=\{\underbrace{\text{LL}}_{\text{approssimazione}},\; \underbrace{\text{LH}}_{\text{dettagli orizzontali}},\; \underbrace{\text{HL}}_{\text{dettagli verticali}},\; \underbrace{\text{HH}}_{\text{dettagli diagonali}}\}. \]
| Sottobanda | Filtri applicati | Contenuto visivo |
|---|---|---|
| LL | basso × basso | Approssimazione dell’immagine (versione smussata e ridotta) |
| LH | basso × alto | Bordi orizzontali e variazioni verticali |
| HL | alto × basso | Bordi verticali e variazioni orizzontali |
| HH | alto × alto | Dettagli diagonali e trame |
La decomposizione può essere applicata ricorsivamente sulla sottobanda LL, generando una rappresentazione multirisoluzione. Dopo \(J\) livelli, si ottiene una struttura con \(3J+1\) sottobande, in cui ogni nuovo livello riduce la risoluzione della componente di approssimazione.
La decomposizione multirisoluzione delle wavelet possiede una relazione concettuale con le rappresentazioni gerarchiche utilizzate nelle reti neurali convoluzionali (CNN). In entrambi i casi, successive fasi di filtraggio e riduzione della risoluzione producono descrizioni sempre più astratte dell’immagine. Tuttavia, le wavelet utilizzano filtri matematicamente definiti e ricostruibili, mentre le CNN apprendono i propri filtri durante l’addestramento.
Diverse famiglie di wavelet presentano compromessi distinti tra supporto spaziale, regolarità e capacità di compressione. Il supporto corrisponde all’estensione della funzione wavelet nel dominio spaziale: quanto più piccolo è il supporto, tanto più localizzata è la funzione; quanto più grande, tanto più liscia tende a essere la sua rappresentazione, sebbene con un costo computazionale maggiore. La Tabella 5.4 confronta alcune delle famiglie più utilizzate.
| Wavelet | Lunghezza del supporto | Momenti nulli | Simmetria | Uso tipico |
|---|---|---|---|---|
| Haar | 2 | 1 | Asimmetrica | Introduzione e analisi di base |
| Daubechies db4 | 8 | 4 | Asimmetrica | Compressione e analisi generale |
| Symlet sym4 | 8 | 4 | Quasi simmetrica | Ricostruzione di segnali |
| Biortogonale 5/3 | 5/3 | 2/2 | Simmetrica | JPEG 2000 senza perdita |
| Biortogonale 9/7 | 9/7 | 4/4 | Simmetrica | JPEG 2000 con perdita |
I momenti nulli misurano la capacità della wavelet di rappresentare regioni uniformi dell’immagine con pochi coefficienti diversi da zero. Una wavelet con \(p\) momenti nulli annulla esattamente i polinomi di grado fino a \(p-1\). Di conseguenza, quanto maggiore è il numero di momenti nulli, tanto maggiore tende a essere l’efficienza di compressione nelle regioni omogenee, sebbene ciò implichi generalmente funzioni con supporto più esteso.
La Figura 5.17 presenta le funzioni di base (wavelet) \(\psi(t)\) nel dominio spaziale. Queste funzioni possiedono supporto compatto, cioè sono diverse da zero solo in una regione finita del dominio, contrariamente alle sinusoidi della Trasformata di Fourier, che si estendono su tutto il dominio.
# ── 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()
Il diagrama della Figura 5.18 illustra l’analisi multirisoluzione eseguita dalla DWT, in cui la sottobanda di approssimazione (LL) viene successivamente decomposta, formando una rappresentazione gerarchica su due livelli.
Il simulatore Figura 5.19 consente di esplorare in modo interattivo la Trasformata Wavelet Discreta 2D (DWT) utilizzando la wavelet di Haar. La decomposizione in sottobande evidenzia la separazione tra la componente di approssimazione e le componenti di dettaglio dell’immagine.
I diversi modelli di ingresso permettono di osservare il comportamento direzionale dei filtri. In immagini con bordi orizzontali e verticali, le sottobande LH e HL evidenziano, rispettivamente, le variazioni verticali e orizzontali dell’intensità. Nelle regioni a variazione graduale, la maggior parte dell’energia si concentra nella sottobanda di approssimazione LL, mentre le sottobande di dettaglio presentano coefficienti prossimi allo zero.
L’analisi multirisoluzione può essere osservata anche aumentando il numero di livelli di decomposizione. In tal caso, solo la sottobanda \(\text{LL}_1\) viene nuovamente decomposta, generando le sottobande \(\text{LL}_2\), \(\text{LH}_2\), \(\text{HL}_2\) e \(\text{HH}_2\), che formano il secondo livello della rappresentazione gerarchica.
In pattern costituiti da regioni omogenee di grande estensione, come una sfumatura graduale o una scacchiera composta da blocchi di grandi dimensioni, l’energia rimane prevalentemente concentrata nella sottobanda LL. Nella sfumatura, ciò avviene perché le differenze tra pixel adiacenti sono piccole. Nella scacchiera, invece, i pixel possiedono praticamente la stessa intensità all’interno di ciascun blocco, cosicché solo i confini tra blocchi producono coefficienti non nulli nelle sottobande di dettaglio. Poiché tali confini occupano solo una piccola frazione dell’immagine, il loro contributo all’energia totale rimane ridotto.
Per rendere possibile l’analisi visiva di queste variazioni sottili, il simulatore incorpora un controllo del guadagno di contrasto dei dettagli (variabile da 1 a 8). Questo parametro funziona come un fattore di amplificazione lineare applicato esclusivamente ai coefficienti delle sottobande di dettaglio (LH, HL e HH) prima della loro visualizzazione a schermo. Negli scenari di transizione graduale (come il gradiente) o di uniformità locale (come l’interno dei blocchi della scacchiera), le differenze numeriche calcolate dal filtro passa-alto di Haar producono coefficienti molto prossimi allo zero, il che renderebbe i quadranti corrispondenti scuri e impercettibili a occhio nudo. Moltiplicando questi valori per il guadagno, il simulatore recupera visivamente le strutture ad alta frequenza nascoste ed enfatizza l’orientamento dei bordi residui.
Il grafico dell’energia per sottobanda quantifica tale distribuzione tra la componente di approssimazione e le componenti di dettaglio, dimostrando che il guadagno visivo non altera la metrica originale dell’energia. Nelle immagini naturali, la maggior parte dell’energia si concentra nella sottobanda LL, mentre le sottobande LH, HL e HH rappresentano principalmente bordi, texture e altre variazioni locali dell’intensità.
La Trasformata Wavelet Discreta 2D (DWT) decompone un’immagine in componenti di approssimazione e dettaglio, organizzate in modo gerarchico su diverse scale e orientazioni. Poiché le sottobande di dettaglio in immagini naturali spesso presentano coefficienti a basso contrasto, gli esempi pratici che seguono utilizzano un pattern geometrico sintetico generato in Python. Questo approccio replica il comportamento del simulatore della Figura 5.19, rendendo visivamente espliciti gli effetti del filtraggio spaziale e della decomposizione multirisoluzione.
La Figura 5.20 illustra la struttura gerarchica della DWT su due livelli utilizzando la wavelet di Haar. Il processo si basa sull’applicazione combinata di filtri passa-basso e passa-alto nelle direzioni orizzontale e verticale, seguiti da un sottocampionamento con fattore 2.
Nel primo livello, l’immagine originale genera la sottobanda di approssimazione (\(LL_1\)) e le componenti di dettaglio orizzontale (\(LH_1\)), verticale (\(HL_1\)) e diagonale (\(HH_1\)). Nell’analisi multirisoluzione, la sottobanda \(LL_1\) viene nuovamente filtrata e sottocampionata, producendo il secondo livello di scomposizione (\(LL_2\), \(LH_2\), \(HL_2\) e \(HH_2\)).
Per consentire l’interpretazione visiva delle componenti di dettaglio, il codice estrae il valore assoluto dei loro coefficienti e applica una normalizzazione lineare (min-max) per occupare l’intera gamma dinamica dei toni di grigio [0, 255]. Questa operazione trasforma le regioni omogenee (coefficienti nulli) in nero ed evidenzia in bianco i bordi e le texture estratte a ciascuna scala e orientazione.
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
# ── Geração da Imagem Sintética (Mesmo padrão 'combined' do simulador) ────────
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)
# Quadrado
if 24 < x < 100 and 24 < y < 100:
v = 225
# Círculo
cx, cy, r = 190, 76, 34
if (x - cx)**2 + (y - cy)**2 < r**2:
v = 205
# Textura periódica (inferior)
if y > 164 and y < 244:
p = 12
v = 185 if ((x // p + y // p) % 2 == 0) else 65
# Linha diagonal
if abs(x - y) < 4:
v = 240
img[y, x] = np.clip(v, 0, 255)
return img.astype(np.uint8)
# Substitui a imagem escura de moedas pelo padrão sintético claro
img_gray = gerar_imagem_sintetica(256)
# ── Decomposição wavelet 2 níveis ─────────────────────────────────────────────
wavelet = "haar"
img_float = img_gray.astype(np.float64)
# Nível 1
coefs1 = pywt.dwt2(img_float, wavelet)
LL1, (LH1, HL1, HH1) = coefs1
# Nível 2 (aplicado sobre LL1)
coefs2 = pywt.dwt2(LL1, wavelet)
LL2, (LH2, HL2, HH2) = coefs2
def sb_vis(sb):
"""Normalizza la sottobanda per la visualizzazione [0,255]."""
return cv2.normalize(np.abs(sb), None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)
print(f"Forma originale : {img_gray.shape}")
print(f"LL1 (livello 1) : {LL1.shape} | LH1/HL1/HH1: {LH1.shape}")
print(f"LL2 (livello 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))Forma originale : (256, 256)
LL1 (livello 1) : (128, 128) | LH1/HL1/HH1: (128, 128)
LL2 (livello 2) : (64, 64) | LH2/HL2/HH2: (64, 64)
La scelta della funzione di base (wavelet) influenza direttamente il modo in cui le caratteristiche dell’immagine vengono distribuite e codificate dai coefficienti della DWT. La Figura 5.21 confronta i risultati pratici ottenuti applicando quattro famiglie distinte sul pattern geometrico sintetico: haar, db4, sym4 e bior2.2.
Poiché possiede un supporto corto e una forma a funzione a gradino, la wavelet di Haar produce coefficienti altamente localizzati in corrispondenza delle discontinuità spaziali, generando bordi sottili e nitidi nelle sottobande di dettaglio. Al contrario, famiglie come Daubechies (db4) e Symlets (sym4), che presentano un supporto maggiore (filtri più lunghi) e un numero più elevato di momenti nulli, generano risposte più morbide e distribuite attorno alle transizioni, il che può introdurre lievi oscillazioni o effetti di sfocatura sui confini netti.
Questo comportamento evidenzia il classico compromesso (trade-off) dell’analisi multirisoluzione: supporti più piccoli favoriscono la localizzazione spaziale esatta dei bordi, mentre supporti più ampi e un maggior numero di momenti nulli tendono a produrre rappresentazioni più sparse e morbide. Questa morbidezza e la capacità di attenuazione delle alte frequenze garantiscono una maggiore efficienza nella compattazione dell’energia, caratteristiche fondamentali per applicazioni di compressione dei dati e rimozione del rumore (denoising).
import numpy as np
import cv2
import pywt
# Garantisce che img_gray e img_float utilizzino lo stesso motivo sintetico chiaro
if 'gerar_imagem_sintetica' in globals():
img_gray = gerar_imagem_sintetica(256)
else:
# Fallback qualora il blocco precedente non sia stato eseguito nella stessa sessione
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))
Una delle principali applicazioni della Trasformata Wavelet Discreta (DWT) è la compressione dei dati, favorita dalla capacità di rappresentazione sparsa dei coefficienti. La Figura 5.22 illustra l’effetto della limiarizzazione netta (hard thresholding), tecnica in cui i coefficienti di dettaglio con magnitudine inferiore a una soglia \(T\) vengono integralmente azzerati prima del processo di sintesi eseguito dalla Trasformata Wavelet Discreta Inversa (IDWT).
Man mano che la soglia \(T\) viene aumentata, un volume crescente di coefficienti ad alta frequenza viene azzerato. Poiché concentrano meno energia, la rimozione di queste componenti riduce considerevolmente la quantità di informazione necessaria per rappresentare l’immagine, mantenendo intatta la componente di approssimazione globale (la sottobanda \(LL\) più profonda) per preservare la struttura macro. Visivamente, questo scarto di coefficienti si manifesta attraverso la scomparsa progressiva delle trame fini e la levigatura delle transizioni brusche di intensità.
La fedeltà dell’immagine ricostruita rispetto all’originale è quantificata dalla metrica del Picco del Rapporto Segnale-Rumore (PSNR, Peak Signal-to-Noise Ratio), espressa in decibel (dB). Valori più elevati di PSNR indicano una distorsione minore e una maggiore prossimità matematica al segnale originale. L’esperimento pratico evidenzia il decadimento graduale del PSNR all’aumentare dell’aggressività della limiarizzazione, consentendo di valutare numericamente la soglia ottimale per il bilanciamento tra compressione e degrado visivo.
import numpy as np
import cv2
import pywt
# Garantisce che img_gray utilizzi lo stesso pattern sintetico chiaro
if 'gerar_imagem_sintetica' in globals():
img_gray = gerar_imagem_sintetica(256)
else:
# Fallback nel caso il blocco precedente non sia stato eseguito nella stessa sessione
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):
"""Decompone, applica la soglia e ricostruisce tramite IDWT."""
coefs = pywt.wavedec2(img.astype(np.float64), wavelet, level=nivel)
# Copia e applica hard thresholding su tutti i dettagli
coefs_t = [coefs[0]] # LL finale non viene sogliato
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)
# Ritaglio alla dimensione originale
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))
La Tabella 5.5 sintetizza le principali differenze strutturali e operative tra la Trasformata Discreta di Fourier (DFT) e la Trasformata Wavelet Discreta (DWT).
| Criterio | Fourier (DFT) | Wavelet (DWT) |
|---|---|---|
| Funzioni di base | Sinusoidi di supporto infinito | Funzioni di supporto compatto |
| Localizzazione spaziale | Non esplicita (globale) | Esplicita (locale) |
| Filtraggio spettrale | Eccellente per il controllo fine delle frequenze | Basata su sottobande (scale) |
| Compressione delle immagini | Base della DCT (JPEG tradizionale) | Base della DWT (JPEG 2000) |
| Analisi multiscala | No | Sì |
| Rimozione del rumore periodico | Altamente efficiente | Poco indicata |
| Segnali non stazionari | Limitata | Altamente efficiente |
In termini pratici, la DFT si consolida come lo strumento ideale per l’analisi spettrale pura, la progettazione di filtri selettivi nel dominio della frequenza e l’attenuazione di rumori periodici e armonici. D’altro canto, la DWT eccelle in scenari che richiedono la rigorosa preservazione della localizzazione spaziale delle caratteristiche associata al loro contenuto frequenziale, distinguendosi nella compressione dei dati, nell’analisi multirisoluzione e nell’elaborazione di transizioni brusche. Pertanto, entrambe le trasformate devono essere comprese come tecniche perfettamente complementari, che tracciano percorsi distinti e specifici per la risoluzione di problemi nell’ambito dell’elaborazione digitale di immagini e visione artificiale (PDI-VC).
Nel tracciare analogie tra l’elaborazione delle immagini e l’audio, è importante considerare le differenze fondamentali:
Nei sistemi audio stereo/multicanale, la fase tra i canali è cruciale per la percezione della localizzazione spaziale (differenze interaurali di fase e di tempo).
Nei sistemi monoaurali, la fase ha un’influenza percettiva limitata — l’orecchio umano è relativamente insensibile alla fase assoluta di componenti sinusoidali isolate.
Nelle immagini, la fase della DFT è sempre fondamentale per la localizzazione spaziale delle strutture, indipendentemente dal fatto che l’immagine sia monocromatica o a colori.
L’analogia tra la fase nell’audio e la fase nelle immagini deve essere utilizzata con cautela, evidenziando che, sebbene entrambe trasportino informazioni sull’organizzazione spaziale/temporale del segnale, i meccanismi percettivi sono fondamentalmente differenti.
Mentre le wavelet stabiliscono il fondamento teorico dello standard JPEG 2000, lo standard JPEG tradizionale si basa sulla Trasformata Discreta dei Coseni (DCT, Discrete Cosine Transform). Nonostante le differenze strutturali, entrambi gli approcci condividono lo stesso principio fondamentale: compattare l’energia dell’immagine in un numero ridotto di coefficienti e scartare le componenti di minore rilevanza con impatto visivo minimo.
L’obiettivo centrale della compressione è ridurre il volume di dati necessario per l’archiviazione o la trasmissione di un’immagine. Tale processo è reso possibile dall’identificazione e dall’eliminazione di ridondanze strutturali e percettive.
Lo sviluppo di algoritmi di compressione si fonda sull’identificazione e sull’eliminazione di tre categorie principali di ridondanza, sintetizzate nella Tabella 5.6.
| Tipo | Definizione | Approccio di Sfruttamento |
|---|---|---|
| Spaziale (interpixel) | Elevata correlazione e dipendenza statistica tra pixel adiacenti. | DCT, DWT e codifica predittiva. |
| Spettrale (intercanale) | Correlazione statistica tra i canali di colore di una stessa immagine. | Trasformazioni dello spazio colore (es: RGB in \(YC_bC_r\)). |
| Psicovisuale | Insensibilità del sistema visivo umano (SVH) alle variazioni ad alta frequenza e a basso contrasto. | Processi di quantizzazione selettiva dei coefficienti. |
A seconda della preservazione dell’informazione originale dopo il processo di decodifica, i metodi di compressione si dividono in due classi fondamentali:
La Trasformata del Coseno Discreta (DCT) costituisce l’operazione centrale dello standard JPEG. A differenza della DFT, che utilizza una base complessa, la DCT si basa su funzioni trigonometriche puramente reali. Per un blocco immagine \(f(x,y)\) di dimensioni \(N \times N\), la DCT-II 2D mappa il segnale spaziale nel dominio delle frequenze spaziali, generando la matrice dei coefficienti \(C(u,v)\) tramite:
\[ 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}\]
dove i fattori di normalizzazione ortogonale sono dati da \(\alpha(0) = \sqrt{1/N}\) e \(\alpha(k) = \sqrt{2/N}\) per \(k > 0\).
Ogni coefficiente \(C(u,v)\) quantifica il contributo — o “peso” — di una specifica frequenza spaziale all’interno di quel blocco. Il termine \(C(0,0)\) è denominato componente DC e rappresenta l’intensità media del blocco (frequenza nulla). I restanti coefficienti, chiamati componenti AC (Alternating Current), corrispondono a frequenze spaziali progressivamente più elevate.
Da una prospettiva geometrica, la Equazione 5.9 realizza la proiezione del blocco di pixel su un insieme di funzioni ortogonali. Per il caso standard del JPEG (\(N=8\)), il blocco spaziale viene decomposto in una combinazione lineare di 64 funzioni di base bidimensionali, denotate da \(B_{u,v}(x,y)\) e generate dal prodotto di funzioni cosinusoidali:
\[B_{u,v}(x,y) = \cos\left[ \frac{\pi (2x+1)u}{16} \right] \cos\left[ \frac{\pi (2y+1)v}{16} \right]\]
In questo modo, l’operazione inversa può essere interpretata come la ricostruzione esatta del blocco originale tramite la somma ponderata di queste 64 matrici di base, dove ogni coefficiente \(C(u,v)\) funge da peso analitico della rispettiva componente armonica.
La frequenza spaziale indicata dagli indici \((u,v)\) determina il numero di cicli di oscillazione lungo le dimensioni orizzontali e verticali del blocco. Come illustrato nella Figura 5.23 — il cui codice isola ciascuna base applicando la trasformazione inversa su impulsi unitari —, queste 64 funzioni sono organizzate in una matrice \(8 \times 8\). L’angolo superiore sinistro (\(u=0, v=0\)) mostra il pattern uniforme a frequenza nulla (DC), mentre il progredire verso destra (asse \(u\)) o verso il basso (asse \(v\)) mappa variazioni armoniche progressivamente maggiori, rappresentando transizioni rapide, bordi e trame nelle orientazioni orizzontali, verticali e diagonali.
Sia la DCT che la DFT mappano un blocco spaziale \(N \times N\) in una matrice di coefficienti della stessa dimensione. Tuttavia, per immagini naturali, la DCT presenta una maggiore efficienza nella compattazione dell’energia alle basse frequenze. Ciò avviene perché la DCT assume implicitamente una simmetria pari del segnale ai bordi del blocco, il che equivale a un’estensione periodica continua, minimizzando l’effetto di dispersione spettrale (ringing). Di conseguenza, la maggior parte dei coefficienti AC decade rapidamente verso valori prossimi allo zero, ottimizzando il pipeline di compressione senza introdurre una degradazione visiva percepibile.
from scipy.fft import dct, idct # riga aggiunta
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()
Prima dell’applicazione della DCT, i pixel del blocco di intensità vengono di routine traslati (sottraendo \(128\) per immagini a 8 bit) al fine di centrare il segnale attorno allo zero, eliminando componenti continue superflue. Calcolando la DCT sul blocco risultante, la proprietà di compattazione dell’energia diventa evidente: la quasi totalità della varianza e dell’informazione dell’immagine originale si concentra nel coefficiente DC (\(C(0,0)\)) e nei primi armonici AC a bassa frequenza.
La Figura 5.24 dimostra questo fenomeno mediante una ricostruzione progressiva per troncamento brusco. Invece di utilizzare tutti i 64 coefficienti, l’algoritmo conserva solo i primi \(k\) componenti — selezionati in base a una scansione che privilegia le basse frequenze spaziali — e azzera i rimanenti.
La sintesi inversa (IDCT) eseguita con solo una frazione dei coefficienti (come il 15% o il 30%) è già in grado di recuperare le strutture e l’illuminazione macro del blocco originale di pixel. Man mano che gli armonici a frequenze più elevate vengono progressivamente reintegrati, i dettagli fini e le transizioni rapide vengono ripristinati. Questo comportamento valida il principio della compressione percettiva: le alte frequenze scartate possiedono poca energia e la loro assenza, in condizioni normali, genera un impatto visivo secondario sulla percezione dell’osservatore.
from scipy.fft import dct, idct
def dct2(bloco):
"""DCT-II 2D ortogonale (separabile)."""
return dct(dct(bloco.T, norm='ortho').T, norm='ortho')
def idct2(coefs):
"""IDCT-II 2D ortogonale."""
return idct(idct(coefs.T, norm='ortho').T, norm='ortho')
# ── Blocco 8×8 centralizzato dell'immagine ─────────────────────────────────────────
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("Coefficienti DCT del blocco 8×8:")
print(np.round(C).astype(int))
print(f"\nEnergia DC : {C[0,0]**2:.1f}")
print(f"Energia totale : {(C**2).sum():.1f}")
print(f"Frazione nel DC : {C[0,0]**2 / (C**2).sum():.1%} ← concentrazione di energia")
# ── Ricostruzione progressiva ──────────────────────────────────────────────────
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))Coefficienti DCT del blocco 8×8:
[[ 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]]
Energia DC : 202725.1
Energia totale : 639868.0
Frazione nel DC : 31.7% ← concentrazione di energia
Lo standard JPEG opera dividendo l’immagine in blocchi disgiunti di \(8 \times 8\) pixel, elaborati tramite una sequenza di trasformazioni spaziali, percettive e statistiche. Il pipeline completo di codifica è strutturato in sei fasi principali:
\[ \text{RGB} \xrightarrow{\text{(1) } YC_bC_r} \xrightarrow{\text{(2) Sottocampionamento}} \xrightarrow{\text{(3) Blocchi } 8 \times 8} \xrightarrow{\text{(4) DCT}} \xrightarrow{\text{(5) Quantizzazione}} \xrightarrow{\text{(6) Codifica Entropica}} \]
La Tabella 5.7 dettaglia la funzione analitica e il fondamento percettivo che giustifica ciascuna di queste fasi.
| Fase | Operazione | Fondamento Percettivo e Statistico |
|---|---|---|
| 1 | Conversione \(RGB \rightarrow YC_bC_r\) | Separa la luminanza (\(Y\)) dalla crominanza (\(C_b, C_r\)). Il sistema visivo umano (SVH) presenta maggiore sensibilità alle variazioni di luminosità che di colore. |
| 2 | Sottocampionamento della crominanza (es: 4:2:0) | Riduce la risoluzione spaziale dei canali colore della metà, scartando dati ridondanti con impatto visivo trascurabile. |
| 3–4 | Centratura e applicazione della DCT \(8 \times 8\) | Trasla i pixel nell’intervallo \([-128, 127]\) e compatta l’energia spettrale del blocco nei coefficienti a bassa frequenza. |
| 5 | Quantizzazione lineare selettiva | Divide ciascun coefficiente \(C(u,v)\) per l’elemento corrispondente della matrice \(Q(u,v)\), applicando un arrotondamento a interi. Costituisce la principale fonte di compressione con perdita. |
| 6 | Scansione a zig-zag e codifica | Ordina i coefficienti quantizzati per massimizzare le sequenze nulle consecutive, ottimizzando la codifica a lunghezza di corsa (RLE) e la codifica di Huffman. |
La matrice di quantizzazione \(Q(u,v)\) è il meccanismo centrale di controllo del compromesso tra tasso di compressione e qualità visiva. Nell’algoritmo pratico della Figura 5.25, il fattore di qualità stabilito dall’utente (scala da 1 a 100) viene convertito in uno scalare che parametrizza la severità della matrice \(Q\). Valori ridotti di qualità espandono i divisori di \(Q(u,v)\), forzando il troncamento di massa dei coefficienti AC a zero. Quando questa eliminazione è eccessiva, la discontinuità ai confini dei blocchi adiacenti non viene attenuata nella ricostruzione, generando i cosiddetti artefatti a blocchi (blocking artifacts).
L’efficienza del codificatore entropico successivo alla quantizzazione dipende direttamente dall’ordinamento dei dati. Poiché la DCT concentra l’energia vitale nel vertice superiore sinistro della matrice (basse frequenze) e spinge i coefficienti nulli verso le estremità opposte, la lettura lineare per righe o colonne frammenterebbe le sequenze di zeri.
L’ordinamento a zig-zag risolve questa limitazione percorrendo la matrice diagonalmente in ordine crescente di frequenza spaziale. Questa mappatura raggruppa i coefficienti significativi all’inizio del vettore e concentra i coefficienti nulli in un’unica sequenza continua alla fine dell’arrangiamento, consentendo all’algoritmo RLE di codificare grandi blocchi di dati in modo compatto ed efficiente.
RLE (Run-Length Encoding) è una tecnica di compressione senza perdita che codifica sequenze consecutive di valori identici — specialmente zeri — come una coppia (conteggio, valore). Nel JPEG, dopo la scansione a zig-zag, i coefficienti quantizzati sono organizzati in modo che gli zeri si concentrino alla fine del vettore. L’RLE comprime quindi questa lunga corsa di zeri con estrema efficienza, ottimizzando l’archiviazione e la trasmissione dell’immagine compressa.
import numpy as np
import cv2
from scipy.fft import dct, idct
# ── Caricamento Sicuro dell'Immagine della Fotocamera (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()
# Ridimensiona leggermente a 256x256 per mantenere lo standard e la velocità dei test precedenti
img_gray = cv2.resize(img_gray, (256, 256))
# ── Tabella di quantizzazione della luminanza (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):
"""DCT-II 2D ortogonale (separabile)."""
return dct(dct(bloco.T, norm='ortho').T, norm='ortho')
def idct2(coefs):
"""IDCT-II 2D ortogonale."""
return idct(idct(coefs.T, norm='ortho').T, norm='ortho')
def jpeg_compress_block(bloco, Q_table):
"""DCT → quantizzazione → dequantizzazione → IDCT in blocchi 8×8."""
C = dct2(bloco.astype(np.float64) - 128)
Cq = np.round(C / Q_table) * Q_table # quantizza e dequantizza
return np.clip(idct2(Cq) + 128, 0, 255)
def jpeg_quality_compress(img, qualidade=50):
"""JPEG semplificato: comprime l'immagine intera per blocchi 8×8."""
if qualidade < 50:
escala = 5000 / qualidade
else:
escala = 200 - 2 * qualidade
# Corretto da 'scala' a '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)
# ── Confronto dei fattori di qualità ───────────────────────────────────
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))
Il simulatore della Figura 5.26 consente di esplorare l’impatto del processo di quantizzazione su un blocco \(8 \times 8\) estratto da un’immagine reale, sintetizzando in tempo reale le seguenti componenti:
La scelta di un formato di memorizzazione digitale incide direttamente sul compromesso tra qualità visiva, dimensione del file e costo computazionale della decodifica. I tre formati di maggiore rilevanza per architetture web e sistemi di calcolo visivo sono JPEG, PNG e WebP.
La Tabella 5.8 sintetizza le proprietà strutturali dei principali formati di immagine rasterizzati.
| Caratteristica | JPEG | PNG | WebP |
|---|---|---|---|
| Compressione | Con perdita | Senza perdita | Con e senza perdita. |
| Trasparenza (canale alfa) | No | Sì | Sì. |
| Supporto animazioni | No | Limitato (APNG) | Sì. |
| Algoritmo di base | DCT + Huffman | DEFLATE (LZ77 + Huffman) | VP8 / VP8L. |
| Ideale per | Fotografia | Grafica, testo e icone | Uso universale in ambiente Web. |
| Inadeguato per | Testo e bordi netti | Immagini fotografiche complesse | Compatibilità legacy. |
Due metriche oggettive sono ampiamente adottate per quantificare la distorsione introdotta dai processi di compressione:
Picco del Rapporto Segnale-Rumore (PSNR, Peak Signal-to-Noise Ratio): \[ \text{PSNR} = 10\,\log_{10}\!\left(\frac{L^2}{\text{MSE}}\right) \quad [\text{dB}] \tag{5.10}\]
dove \(L = 255\) per immagini quantizzate a 8 bit e \(\text{MSE}\) rappresenta l’Errore Quadratico Medio (Mean Squared Error). Valori di PSNR superiori a 40 dB indicano fedeltà eccellente; tra 30 dB e 40 dB rappresentano buona qualità; valori inferiori a 30 dB corrispondono a degradazioni visive facilmente percepibili.
Indice di Similarità Strutturale (SSIM, Structural Similarity Index): \[ \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}\]
Lo SSIM valuta finestre locali dell’immagine basandosi su tre componenti complementari: luminanza (\(\mu_f, \mu_g\)), contrasto (\(\sigma_f, \sigma_g\)) e struttura (\(\sigma_{fg}\)), ponderate da costanti di stabilità \(c_1\) e \(c_2\). L’indice varia nell’intervallo \([-1, 1]\), dove l’unità rappresenta l’identità perfetta. A differenza del PSNR, lo SSIM considera l’organizzazione spaziale degli errori, allineandosi alla percezione del sistema visivo umano (SVH).
Il PSNR possiede una formulazione matematica semplice e un basso costo computazionale; tuttavia tende a sovrastimare la qualità in immagini con distorsioni localizzate o a sottostimarla in variazioni globali di luminosità tollerate dall’osservatore. Lo SSIM modella con maggiore fedeltà la percezione biologica, ma richiede un maggiore sforzo di elaborazione. Per analisi rigorose dei codec, si raccomanda di riportare entrambe le metriche statistiche in carattere complementare.
La natura matematica del codificatore determina il tipo di degrado introdotto a bitrate ridotti. Come illustrato nella Figura 5.27, la compressione aggressiva tramite DCT nello standard JPEG segmenta l’immagine in griglie rigide, generando gli artefatti a blocchi (blocking artifacts). Al contrario, gli algoritmi basati su codifica predittiva o rappresentazioni sottoposte a trasformate spaziali avanzate (come WebP e JPEG 2000) eliminano le discontinuità di blocco, ma introducono perdita di texture fine e sfocature caratteristiche attorno ai bordi ad alto contrasto.
import os
import cv2
# Garantisce l'esistenza della directory e salva i file compressi
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])
# Estrazione della regione di interesse per la visualizzazione degli artefatti (Zoom 4x)
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=["Zoom Originale", "JPEG Q=10 (Artefatto di Blocco)", "WebP Q=10 (Smussamento)"],
cols=3, figsize=(14, 5))
La validazione degli algoritmi di compressione con perdita richiede un’analisi che correli il costo di archiviazione alla fedeltà del segnale ricostruito. Tale valutazione viene condotta in modo complementare attraverso curve di prestazione globale e mediante la mappatura locale delle distorsioni indotte dai codificatori.
La Figura 5.28 presenta la valutazione empirica del pipeline JPEG e WebP tramite curve di rateo-distorsione, che monitorano il guadagno di compressione (dimensione del file in KB) in funzione del PSNR. Il formato PNG funge da linea di base ideale (\(\text{PSNR} = \infty\)), poiché la sua natura lossless impedisce qualsiasi degradazione, sebbene richieda un volume di dati notevolmente maggiore.
L’analisi delle curve dimostra la superiorità e l’efficienza dello standard WebP rispetto al JPEG tradizionale: per raggiungere lo stesso livello di fedeltà matematica (come la fascia di qualità eccellente, dove \(\text{PSNR} > 40\text{ dB}\)), il codificatore WebP genera file significativamente più piccoli. Questo comportamento riflette l’impatto pratico dell’evoluzione degli algoritmi nell’ottimizzazione dei sistemi di trasmissione e archiviazione digitale.
import os
import cv2
import matplotlib.pyplot as plt
# Garantisce l'esistenza della directory di test
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})
# ── Generazione della Curva Rateo-Distorsione ─────────────────────────────────
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"\nDimensione grezza (senza compressione): {img_gray.nbytes/1024:.0f} KB")
print(f"\n{'Formato':>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}")
Dimensione grezza (senza compressione): 64 KB
Formato 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.14L’immagine Cameraman (\(256 \times 256\) pixel in scala di grigi) occupa 64 KB in formato grezzo (senza compressione). Come riferimento, il PNG lossless comprime questo volume a 36,2 KB — evidenziando che la compressione senza perdita riduce già significativamente lo spazio di archiviazione per immagini con regioni omogenee. In contrapposizione, i formati con perdita (JPEG e WebP) raggiungono dimensioni ancora minori: il JPEG con qualità 95 occupa 22,3 KB (PSNR ≈ 45 dB), mentre il WebP con qualità 90 raggiunge 12,5 KB con PSNR equivalente, dimostrando la sua superiorità in efficienza di compressione.
Sebbene il PSNR offra un indicatore numerico rapido, le metriche globali non riescono a distinguere come la perdita di informazione si distribuisca geometricamente sull’immagine. La Figura 5.29 risolve questa limitazione associando le ricostruzioni a diverse qualità ai rispettivi mappe di errore assoluto e all’SSIM.
Le mappe residue — ottenute dalla differenza assoluta normalizzata tra l’immagine originale e quella compressa — rivelano la firma spaziale intrinseca di ciascuna architettura di codifica:
L’SSIM cattura questa degradazione morfologica in modo molto più sensibile rispetto al PSNR, penalizzando il punteggio finale man mano che l’organizzazione strutturale e le trame fini — alle quali il sistema visivo umano è altamente reattivo — vengono eliminate dal codificatore.
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
# Garantisce l'esistenza della directory di test
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"
# REGISTRAZIONE FORZATA: Genera e registra il JPEG con la qualità attuale nel percorso corretto
img_compactada = jpeg_quality_compress(img_gray, qualidade=q)
cv2.imwrite(path, img_compactada)
# Lettura sicura del file appena registrato
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)
# Differenza assoluta normalizzata per evidenziare la struttura spaziale dell'errore
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))
Le mappe di errore presentate sono state normalizzate singolarmente (cv2.NORM_MINMAX) per massimizzare il contrasto visivo e rivelare la struttura spaziale delle distorsioni. Ciò significa che:
Pertanto, l’intensità del bianco nelle mappe NON è confrontabile tra diverse qualità — le mappe servono solo a rivelare la firma spaziale dell’errore (bordi vs blocchi), non la sua ampiezza. L’ampiezza corretta è fornita dai valori di PSNR e SSIM, che mostrano chiaramente che Q=95 ha un errore molto minore rispetto a Q=25.
Il processo di compressione nello standard JPEG si basa sull’applicazione combinata di trasformazioni spaziali, percettive e statistiche per ridurre le ridondanze di un’immagine. La Tabella 5.9 riassume il ruolo di ciascuna fase nel pipeline e il rispettivo impatto sulla riduzione dei dati.
| Fase | Operazione Analitica | Meccanismo di Guadagno / Compressione |
|---|---|---|
| Conversione \(YC_bC_r\) | Isolamento dei canali di luminanza e crominanza. | Modella la percezione del SVH, consentendo di trattare colore e luminosità in modo indipendente. |
| Sottocampionamento 4:2:0 | Riduzione della risoluzione spaziale dei canali di colore (\(C_b\) e \(C_r\)). | Elimina circa il 50% dei dati grezzi con un impatto visivo minimo. |
| DCT \(8 \times 8\) | Mappatura dal dominio spaziale al dominio delle frequenze spaziali. | Compattazione dell’energia, concentrando l’informazione vitale nei primi coefficienti. |
| Quantizzazione Lineare | Divisione intera dei coefficienti per una matrice di ponderazione \(Q(u,v)\). | Principale fonte di compressione con perdita; elimina le alte frequenze impercettibili. |
| Codifica Entropica | Applicazione di algoritmi RLE e codifica di Huffman. | Compressione statistica senza perdita, ottimizzata dalle lunghe sequenze di coefficienti nulli. |
L’applicazione di tassi di compressione eccessivamente aggressivi (fattori di qualità ridotti) introduce distorsioni prevedibili nell’immagine ricostruita, derivanti dalle limitazioni matematiche del modello:
Integrando le tecniche consolidate nel corso di questo capitolo, si presenta un pipeline completo di restauro delle immagini che combina l’analisi spettrale nel dominio della frequenza con il filtraggio adattivo nel dominio spaziale. L’obiettivo è attenuare un rumore misto (composto da degradazione gaussiana e interferenza periodica) preservando al massimo i dettagli strutturali dell’immagine originale.
\[ \text{Immagine Rumorosa} \xrightarrow{\text{FFT2}} \xrightarrow{\text{Filtro Notch Gaussiano}} \xrightarrow{\text{IFFT2}} \xrightarrow{\text{Filtro Bilaterale}} \text{Immagine Restaurata} \]
La coppia di metriche statistiche PSNR e SSIM fornisce una valutazione qualitativa e morfologica complementare del processo di restauro:
Nella pratica, esiste un compromesso analitico (trade-off) tra riduzione del rumore e preservazione dei dettagli: filtri spaziali eccessivamente aggressivi attenuano bene il rumore ad alta frequenza, ma degradano trame fini e smussano bordi netti — il che riduce simultaneamente sia il PSNR che lo SSIM rispetto all’immagine originale. La sfida nella progettazione dei filtri è trovare il punto di equilibrio che massimizzi entrambe le metriche, garantendo un restauro fedele e visivamente gradevole.
I risultati numerici e visivi generati dalla Figura 5.30 dimostrano la rilevanza pratica di associare diversi domini di elaborazione. L’inserimento simultaneo di rumore periodico e stocastico corrompe le proprietà morfologiche del segnale, riducendo severamente gli indici di similarità e il rapporto segnale-rumore dell’immagine di riferimento.
L’isolamento e la soppressione dei picchi armonici nel dominio della frequenza tramite la maschera notch rimuovono le frange di interferenza sinusoidali sparse nello spazio bidimensionale. Come evidenziato nei dati stampati della Figura 5.30, questa filtrazione chirurgica promuove un salto immediato e sostanziale nella metrica PSNR. Tuttavia, il rumore gaussiano ad alta frequenza rimane attivo in modo omogeneo nello spettro, richiedendo un approccio complementare.
Il restauro finale è consolidato nel dominio spaziale con l’introduzione del filtro bilaterale. Diversamente dagli operatori passa-basso convenzionali (come quello gaussiano o di media), che smusserebbero indiscriminatamente rumore e contorni strutturali, la filtrazione bilaterale calcola pesi ponderati in base alla prossimità geometrica e alla differenza di intensità radiometrica. Questo comportamento adattivo attenua le fluttuazioni stocastiche residue nelle regioni di transizione graduale e preserva la nitidezza dei bordi spaziali.
La convergenza di entrambi gli approcci risulta in un miglioramento sostanziale e simultaneo di PSNR e SSIM rispetto all’immagine rumorosa — sebbene i valori finali rimangano inferiori a quelli dell’immagine originale (PSNR = \(\infty\), SSIM = 1,0), a causa della perdita inevitabile di informazioni spettrali e testurali durante i processi di filtrazione. L’attenuazione graduale (gaussiana) dei picchi nello spettro evita artefatti di ringing, mentre il filtro bilaterale elimina il rumore stocastico residuo senza compromettere la nitidezza dei bordi. I risultati confermano l’efficacia e la complementarità pratica degli strumenti di analisi di frequenza presentati in questo capitolo, dimostrando che la filtrazione ibrida (frequenza + spaziale) è superiore a qualsiasi approccio isolato per il restauro di immagini degradate da rumore misto.
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. Costruzione del rumore misto ─────────────────────────────────────────────
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. Spettro e identificazione dei picchi ────────────────────────────────────
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. Notch gaussiano sui picchi periodici ──────────────────────────────────
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. Filtro bilaterale: rimozione del rumore gaussiano residuo ─────────────────
img_den = cv2.bilateralFilter(img_notch, d=7, sigmaColor=25, sigmaSpace=7)
# ── Calcolo delle metriche di validazione ────────────────────────────────────────
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=[
"Originale",
f"Rumosa\nPSNR={psnr_n:.1f} dB",
"Spettro\n(picchi visibili)",
"Maschera notch\n(gaussiana morbida)",
f"Dopo notch\nPSNR={psnr_no:.1f} dB",
f"Notch + bilaterale\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
La transizione dal dominio spaziale al dominio delle frequenze rivela la distribuzione spettrale dell’energia dell’immagine, stabilendo la base analitica per il filtraggio avanzato, il restauro e la compressione dei dati. L’articolazione strutturale di questi concetti è sintetizzata nella mappa concettuale della Figura 5.31.
Prossimi Passi: Il Capitolo 6 inaugura la Parte II dell’opera, applicando gli strumenti di elaborazione delle immagini alla risoluzione di problemi reali di ispezione industriale. Verranno esplorate tecniche di segmentazione e analisi delle forme per il rilevamento automatico di difetti nelle linee di produzione — dall’identificazione di difetti superficiali nei pezzi alla lettura di QRCode nelle prove, consolidando il ponte tra la teoria presentata nella Parte I e le esigenze pratiche della visione computazionale.
In questa edizione, si incoraggia l’uso della piattaforma Gemini Notebook come strumento complementare di apprendimento — non come sostituto della lettura attenta, della risoluzione degli esercizi o della sperimentazione pratica. Basato su architetture di intelligenza artificiale, il sistema utilizza esclusivamente il materiale didattico e i documenti forniti dall’autore come base di conoscenza, garantendo che le risposte generate siano concettualmente allineate al contenuto programmatico e all’approccio pedagogico adottato nel corso di quest’opera.
(10%) Implementazione Diretta della DFT 2D: Implementare analiticamente la Trasformata Discreta di Fourier 2D (DFT) senza l’ausilio di funzioni native di librerie (come np.fft.fft2), utilizzando strettamente la formulazione matematica definita in Equazione 5.1 per una matrice di dimensioni \(16 \times 16\). Eseguire la validazione numerica confrontando i coefficienti generati con i risultati della funzione np.fft.fft2, assicurandosi che la deviazione assoluta massima sia inferiore a \(10^{-8}\). Misurare i tempi di esecuzione di entrambi i metodi e presentare una giustificazione teorica per la disparità osservata in termini di complessità asintotica.
(15%) Soppressione del Rumore Periodico: Aggiungere interferenze sinusoidali con frequenze spaziali \((u_0, v_0) \in \{(5,10), (20,5), (30,30)\}\) all’immagine di test del Cameraman. Per ogni scenario di degradazione, progettare una maschera di filtraggio notch specifica nel dominio della frequenza per isolare e attenuare i picchi armonici indesiderati. Valutare quantitativamente l’efficacia del processo di restauro mediante il calcolo delle metriche PSNR e SSIM. Discutere analiticamente il compromesso (trade-off) tra l’attenuazione del rumore sinusoidale e l’indesiderata attenuazione delle caratteristiche strutturali legittime dell’immagine.
(15%) Analisi Comparativa degli Operatori Passa-Basso: Condurre uno studio comparativo tra i filtri passa-basso Ideale, Gaussiano e Butterworth (con ordini armonici \(n = 1, 2, 4\)), parametrizzati con frequenze di taglio \(D_0 = 20, 40, 60\) pixel. Per ogni combinazione strutturale, calcolare gli indici PSNR e SSIM dell’immagine risultante rispetto al segnale originale di riferimento. Organizzare i dati quantitativi in una tabella strutturata e tracciare i grafici unidimensionali delle funzioni di trasferimento corrispondenti lungo il profilo orizzontale \(H(u, 0)\).
(15%) Banco di Filtri Multirisoluzione di Haar: Sviluppare uno script per eseguire manualmente la decomposizione wavelet discreta 2D di primo livello utilizzando la famiglia Haar. L’algoritmo deve calcolare i coefficienti dei filtri corrispondenti passa-basso (\(h\)) e passa-alto (\(g\)), applicandoli in modo separabile sulle righe e colonne della matrice, seguiti dall’operazione di decimazione (sottocampionamento spaziale per un fattore di 2). Validare numericamente l’accuratezza della propria implementazione confrontando le sottobande ottenute con l’output della funzione pywt.dwt2(img, 'haar').
(15%) Compressione Sparsa mediante Sogliatura Wavelet: Applicare la tecnica di filtraggio per sogliatura netta (hard thresholding) sui coefficienti di dettaglio della decomposizione wavelet, adottando le soglie numeriche \(T \in \{5, 10, 20, 40, 80\}\) per le famiglie Haar, Daubechies (db4) e Symlets (sym4). Dopo aver eseguito il processo di sintesi mediante la trasformata inversa (pywt.waverec2), calcolare i valori di PSNR e SSIM di ciascuna immagine ricostruita. Identificare e giustificare quale combinazione di famiglia wavelet e soglia \(T\) massimizza la similarità strutturale.
(15%) Costruzione di un Codificatore JPEG Semplificato: Implementare la pipeline completa di compressione dei dati simulando lo standard JPEG. Il flusso deve comprendere: conversione spaziale \(RGB \rightarrow YC_bC_r\), sottocampionamento cromatico nella proporzione 4:2:0, segmentazione della luminanza in blocchi disgiunti di \(8 \times 8\) pixel, applicazione della DCT-II 2D ortogonale e quantizzazione lineare basata sulla matrice normalizzata di luminanza scalata per i fattori di qualità desiderati. Eseguire la decodifica inversa e confrontare quantitativamente le ricostruzioni con i file generati dalla funzione cv2.imencode per i fattori di qualità di 20, 50 e 80.
(15%) Analisi Percettiva su Contenuti Eterogenei: Sviluppare un’immagine sintetica composta da tre regioni distinte e di caratteristiche spettrali contrastanti: una texture fotografica complessa (che rappresenta alte frequenze stocastiche), un’area di testo vettorizzato con bordi netti (che rappresenta transizioni a gradino pure) e un gradiente lineare continuo (che rappresenta basse frequenze omogenee). Sottoporre questa immagine mista ai processi di compressione nei formati JPEG, PNG e WebP. Valutare e interpretare i risultati correlando la dimensione finale del file su disco alle metriche PSNR e SSIM ottenute, giustificando quale formato mostra le prestazioni migliori per segnali di natura eterogenea e perché tale vantaggio si verifica in termini di compattazione dell’energia e preservazione percettiva.
La base teorica e lo sviluppo analitico dei concetti trattati in questo capitolo si fondano sulle seguenti opere di riferimento:
La presente lista di esercizi di programmazione (EP) consolida le formulazioni teoriche presentate nel corso del Capitolo 5 — Trasformate e Compressione — attraverso un percorso pratico applicato. Gli esercizi sono strutturati a partire da matrici di dimensioni ridotte, consentendo la validazione analitica e l’ispezione manuale di ciascun coefficiente, mantenendo la coerenza metodologica adottata nei capitoli precedenti.
L’incatenamento degli esercizi riproduce rigorosamente il flusso concettuale del capitolo: si inizia con l’implementazione esplicita della Trasformata Discreta di Fourier (DFT) a partire dalla sua definizione matematica fondamentale; si prosegue con la progettazione di filtri passa-basso e maschere notch nel dominio della frequenza; si applica la quantizzazione dei coefficienti (nucleo della compressione con perdita); e si conclude con l’integrazione di queste fasi nella costruzione di un pipeline di compressione JPEG semplificato e nell’analisi percettiva dei formati immagine.
In tutti gli esercizi di questo capitolo, le coordinate del centro dello spettro (origine delle frequenze spaziali dopo l’applicazione dello spostamento fftshift) devono essere determinate tramite divisione intera. Per una matrice con \(L\) righe e \(C\) colonne, la componente di frequenza nulla si trova nella posizione:
\[ (c_y, c_x) = \left( \left\lfloor \frac{L}{2} \right\rfloor, \left\lfloor \frac{C}{2} \right\rfloor \right) \]
Questa convenzione è rigorosamente identica a quella adottata dalla funzione np.fft.fftshift. Inoltre, in tutte le fasi che richiedono discretizzazione o arrotondamento numerico (sia nella quantizzazione dei coefficienti AC sia nella ricostruzione finale dei pixel), si deve impiegare l’arrotondamento standard al numero intero più vicino (round half away from zero), mitigando ambiguità in valori con frazione esattamente pari a \(0.5\).
Il quaderno consente di sviluppare, validare, organizzare e testare soluzioni di Esercizi di Programmazione (EPs) in ambienti interattivi, come Colab, con gli stessi casi di test di Moodle, copiandoli lì solo al momento di registrare il voto ufficiale.
Scarica morph.py e testsuite.py eseguendo la cella qui sotto:
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✅ Ambiente pronto. Morph: 1.1.9 | OpenCV: 5.0.0 | TestSuite: 1.1.2Per valutare i test, esegui TestSuite("EP05_01.estensione").run() in una nuova cella, sostituendo l’estensione con quella del linguaggio utilizzato (.py, .java, .c, .cpp, .js o .r). Il sistema scarica i casi di test da GitHub, esegue il programma e calcola automaticamente il voto.
Per testare il codice Python direttamente, senza salvare un file, usa run_code(codice) passando il codice come stringa in una variabile codice:
codice = """
from morph import mm
# ... il tuo codice qui ...
"""
TestSuite("EP05_01").run_code(codice)In uno scanner di documenti antico, il sensore cattura carta stropicciata e la trama delle fibre insieme al testo — rumore ad alta frequenza che “inquina” lo spettro ai bordi. Il tecnico della manutenzione non ha accesso all’immagine originale, ma solo allo spettro di magnitudo già calcolato dal software dello scanner. Il suo compito è semplice e chirurgico: mantenere solo il cerchio centrale delle basse frequenze (la struttura globale del documento) ed eliminare tutto ciò che si trova al di fuori del raggio \(D_0\), rimuovendo la trama fine senza nemmeno dover toccare l’immagine spaziale.
Questo è il Filtro Passa-Basso Ideale (LPFI): l’operazione spettrale più diretta del capitolo, ma anche quella che meglio rivela l’anatomia di uno spettro centrato.
np.fft.fftshift).| Regione | Distanza dal centro | Effetto del filtro |
|---|---|---|
| Centro (\(D \le D_0\)) | Basse frequenze | Preservate — struttura globale mantenuta |
| Bordi (\(D > D_0\)) | Alte frequenze | Azzerate — trama e rumore rimossi |
| \(D_0\) piccolo | — | L’immagine ricostruita sarebbe molto sfocata |
| \(D_0\) grande | — | Poca filtrazione; quasi tutta l’energia preservata |
Ingresso:
Uscita:
| Ingresso | Uscita | Osservazione |
|---|---|---|
| 3 3 1 10 20 30 40 50 60 70 80 90 |
0 20 0 40 50 60 0 80 0 |
Centro \((1,1)\). Gli angoli hanno \(D=\sqrt{2}\approx1.41 > 1\), quindi vengono azzerati; i vicini ortogonali hanno \(D=1 \le 1\) e sono mantenuti. |
| 1 3 0 5 9 7 |
0 9 0 | \(L=1, C=3\): centro in \((0,1)\). Solo la posizione centrale stessa (\(D=0\)) sopravvive a \(D_0=0\). |
%%writefile EP05_01.py
# Codice PythonOverwriting EP05_01.pyTestSuite("EP05_01.py").run()✔️ EP05_01.cases esiste già in casos/
📋 5 caso/i caricato/i da casos/EP05_01.cases
🔍 Test di Python: EP05_01.py
⚠️ EP05_01.py: file vuoto (meno di 3 righe). Test saltati.
Una telecamera di ispezione industriale acquisisce immagini di circuiti stampati, ma l’alimentazione della linea di produzione introduce un’interferenza elettrica periodica — un pattern di strisce quasi impercettibile a occhio nudo, che però appare nello spettro di Fourier come coppie di picchi luminosi posizionati simmetricamente attorno al centro. Il team di visione artificiale non può rielaborare l’acquisizione: deve localizzare e cancellare chirurgicamente queste coppie di picchi nello spettro, preservando tutta l’altra informazione utile dell’immagine.
Questo è il ruolo del filtro rigetta-banda notch: a differenza del passa-basso (che interessa una regione continua), esso agisce su punti specifici e sui loro simmetrici, lasciando intatto il resto dello spettro.
| Concetto | Ruolo nel filtro notch |
|---|---|
| Picco in \((\Delta v, \Delta u)\) | Frequenza dell’interferenza periodica rilevata visivamente nello spettro |
| Punto simmetrico \((-\Delta v,-\Delta u)\) | Ogni DFT di segnale reale è hermitiana: i picchi compaiono sempre in coppie simmetriche rispetto al centro |
| Raggio \(r\) | Controlla la “larghezza” della reiezione — un \(r\) grande rimuove più energia attorno al picco, ma anche informazione utile |
Input:
Output:
| Input | Output | Osservazione |
|---|---|---|
| 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 |
Centro \((c_y, c_x) = (2, 2)\). Il picco indicato \((\Delta v, \Delta u) = (1, 1)\) genera il punto \((3, 3)\) (valore 19) e il suo simmetrico \((1, 1)\) (valore 7), entrambi azzerati con \(r=0\) (solo i punti esatti). |
%%writefile EP05_02.py
# Codice PythonOverwriting EP05_02.pyTestSuite("EP05_02.py").run()✔️ EP05_02.cases esiste già in casos/
📋 5 caso/i caricato/i da casos/EP05_02.cases
🔍 Test di Python: EP05_02.py
⚠️ EP05_02.py: file vuoto (meno di 3 righe). Test saltati.
Un’applicazione di galleria fotografica deve ridurre le dimensioni di migliaia di immagini prima di caricarle sul cloud, senza ricodificare tutto da zero. L’ingegnere responsabile ha già i coefficienti DCT di ogni blocco \(4\times4\) calcolati (la fase computazionalmente onerosa è già stata eseguita) — manca solo applicare la tabella di quantizzazione, la fase che scarta realmente informazioni e genera compressione. I coefficienti ad alta frequenza, meno percettibili all’occhio umano, ricevono divisori grandi e tendono a diventare zero; i coefficienti a bassa frequenza, più percettibili, ricevono divisori piccoli e sopravvivono quasi intatti.
Dovrai implementare esattamente questa fase: quantizzare e dequantizzare (dividere, arrotondare, moltiplicare di nuovo) — il cuore della compressione lossy del JPEG.
.5 non si verificano mai nei casi di test).| Coefficiente | Frequenza | Valore tipico di \(Q\) | Effetto della quantizzazione |
|---|---|---|---|
| \(C(0,0)\) | DC (media del blocco) | Piccolo | Quasi sempre sopravvive — domina l’energia |
| \(C(u,v)\) con \(u+v\) basso | Bassa frequenza | Piccolo/medio | Parzialmente preservato |
| \(C(u,v)\) con \(u+v\) alto | Alta frequenza | Grande | Spesso diventa zero — fonte della compressione |
Input:
Output:
| Input | Output | Osservazione |
|---|---|---|
| 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\) (preservato). \(C(0,2)=-5/7\approx-0.71\to-1\to-1\times7=-7\). Invece \(C(1,1)=-3/7\approx-0.43\to0\): azzerato dalla quantizzazione — la maggior parte del blocco diventa zero, illustrando la compattazione dell’energia nell’angolo superiore sinistro. |
%%writefile EP05_03.py
# Codice PythonOverwriting EP05_03.pyTestSuite("EP05_03.py").run()✔️ EP05_03.cases esiste già in casos/
📋 5 caso/i caricato/i da casos/EP05_03.cases
🔍 Test di Python: EP05_03.py
⚠️ EP05_03.py: file vuoto (meno di 3 righe). Test saltati.
Un laboratorio di ricerca in astronomia computazionale ha ricevuto, da una missione remota, un piccolo sensore sperimentale i cui dati grezzi non possono essere elaborati tramite librerie moderne di FFT — l’ambiente di validazione è isolato e consente solo operazioni aritmetiche di base. Il team deve reimplementare la Trasformata di Fourier Discreta 2D a partire dalla definizione matematica stessa, cella per cella, per poi confrontare bit per bit con np.fft.fft2 in un altro ambiente.
Questo è l’esercizio più concettuale della lista: non ci sono scorciatoie. Dovrai implementare direttamente la doppia sommatoria della Equazione 5.1, evidenziando perché la FFT esiste — e il costo computazionale che essa evita.
fftshift — il componente DC rimane in \((0,0)\)).fftshift: l’output mantiene la convenzione grezza della DFT, con il componente DC in \(F(0,0)\) (angolo superiore sinistro)..5.| Elemento | Significato |
|---|---|
| \(F(0,0)\) | Componente DC — somma di tutti i pixel, \(F(0,0) = \sum f(x,y)\) |
| Parte reale \(\text{Re}(F)\) | Proiezione del segnale sui coseni |
| Parte immaginaria \(\text{Im}(F)\) | Proiezione del segnale sui seni |
| Complessità di questa implementazione | \(\mathcal{O}((MN)^2)\) — ecco perché la FFT, con \(\mathcal{O}(MN\log(MN))\), è indispensabile nelle immagini reali |
Input:
Output:
| Input | Output | Osservazione |
|---|---|---|
| 2 2 1 2 3 4 |
10 2 4 0 |
\(F(0,0)=1+2+3+4=10\) (DC = somma totale). \(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
# Codice PythonOverwriting EP05_04.pyTestSuite("EP05_04.py").run()✔️ EP05_04.cases esiste già in casos/
📋 5 caso/i caricato/i da casos/EP05_04.cases
🔍 Test di Python: EP05_04.py
⚠️ EP05_04.py: file vuoto (meno di 3 righe). Test saltati.
Sei stato incaricato di creare, da zero, un codec JPEG didattico in un ambiente embedded, senza alcuna libreria di immagini disponibile — solo operazioni matematiche di base. Il cliente vuole capire esattamente dove la qualità viene persa e dove viene recuperata, blocco per blocco. Questa è la sfida finale del capitolo: integrare tutto ciò che è stato studiato — la DCT-II ortonormale, la quantizzazione percettiva e la ricostruzione tramite IDCT — in un unico pipeline end-to-end, elaborando un blocco \(N \times N\) dall’inizio alla fine, esattamente come fa internamente lo standard JPEG, \(8\times8\) pixel alla volta.
.5.| Fase | Analoga nel vero standard JPEG | Dove la qualità viene persa |
|---|---|---|
| Centratura | Stessa — la DCT presuppone un segnale centrato su zero | Nessuna perdita |
| DCT-II | Fasi 3–4 del pipeline (Tabella 5.7) | Nessuna perdita (trasformazione esatta e reversibile) |
| Quantizzazione | Fase 5 — divisione per \(Q(u,v)\) | Principale fonte di perdita — i coefficienti ad alta frequenza diventano zero |
| IDCT | Ricostruzione finale | Ricostruisce esattamente i coefficienti quantizzati, non quelli originali |
Input:
Output:
| Input | Output | Osservazione |
|---|---|---|
| 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 |
Dopo la DCT, una quantizzazione aggressiva sulle alte frequenze (valori grandi di \(Q\) nell’angolo in basso a destra) e la ricostruzione tramite IDCT, il blocco risulta vicino all’originale, ma non identico — la differenza è il costo della compressione lossy. |
Se il risultato non corrisponde, verifica in questo ordine: (1) i coefficienti DCT grezzi (prima della quantizzazione) — devono ricostruire l’originale esattamente tramite IDCT se salti le fasi 6–7; (2) la tabella \(\alpha(u)\) — un errore comune è applicare \(\sqrt{2/N}\) anche per \(u=0\); (3) l’arrotondamento della quantizzazione, che deve avvenire prima di moltiplicare di nuovo per \(Q\).
%%writefile EP05_05.py
# Codice PythonOverwriting EP05_05.pyTestSuite("EP05_05.py").run()✔️ EP05_05.cases esiste già in casos/
📋 5 caso/i caricato/i da casos/EP05_05.cases
🔍 Test di Python: EP05_05.py
⚠️ EP05_05.py: file vuoto (meno di 3 righe). Test saltati.