Esercizio - Laboratorio 4 - interpolazione e zoom di immagini

Testo (laboratorio 4 del corso Multimedia Signal Processing, UniPD, lezione 21; con codice Python numpy/scipy). Il laboratorio applica il filtraggio FIR all'ingrandimento (zooming) di immagini, dove un filtro passabasso fa l'interpolazione:

1. Riscaldamento (1D): generare una sequenza di campioni casuali con due zeri fra un campione non nullo e il successivo ($L=3$), interpolarla linearmente con un FIR a coefficienti triangolari e confrontare il risultato con l'originale: c'è una traslazione? di quanti campioni?
2. Zoom di un'immagine di $50\times50$ pixel di un fattore $5$: (a) ripetizione dei pixel; (b) inserzione di zeri nelle righe e nelle colonne, poi (c) interpolazione lineare con filtro triangolare ($L=5$) separabile; (d) interpolazione con un'approssimazione del filtro passabasso ideale, il sinc campionato $b_k=\operatorname{sinc}(k/5)$, $k=-10,\dots,10$, con e senza finestra di Hamming.
3. Calcolare e rappresentare la risposta in frequenza dei filtri 1D usati.

Dati. L'immagine è Lena.tiff (nella cartella del laboratorio 4 su Moodle; la consegna la chiama Lena.tif): $512\times512$ pixel a colori, che si converte in grigi e in valori $[0,1]$ (livello medio $0{,}486$, minimo $0{,}098$, massimo $0{,}961$). La patch della consegna è img[240:290, 240:290], una finestra di $50\times50$ pixel che contiene un occhio. Teoria usata: Elaborazione multirate - decimazione e interpolazioneL'elaborazione multirate cambia la frequenza di campionamento di un segnale. Interpolazione per $L$ (da $F_s$ a $LF_s$): l'espansore inserisce $L-1$ zeri fra i campioni, $W(\hat\omega)=X(L\hat\omega)$, lo spettro si ripete $L$ volte nell'asse più ampio; un passabasso con taglio $\pi/L$ e guadagno $L$ elimina le immagini, e con il sinc ideale $h[n]=\operatorname{sinc}(n/L)$ i campioni originali sono conservati ($y[nL]=x[n]$); in pratica FIR a finestre o interpolatore lineare (triangolo). Decimazione per $M$ (da $F_s$ a $F_s/M$): $y[n]=x[nM]$, $Y(\hat\omega)=\frac1M\sum_{k=0}^{M-1}X\big(\frac{\hat\omega-2\pi k}M\big)$, repliche sovrapposte (aliasing) se la banda supera $\pi/M$: si fa precedere da un passabasso con taglio $\pi/M$. Conversione razionale $F_s'=\frac LMF_s$ ($L,M$ coprimi): espansore $L$, un solo passabasso di taglio $\min(\frac\pi L,\frac\pi M)$ e guadagno $L$, decimatore $M$; $y[n]=\sum_kx[k],g[nM-kL]$ non è una convoluzione e il sistema è periodicamente tempo-invariante. Realizzare il tutto in più stadi (fattori piccoli prima per l'interpolazione, grandi prima per la decimazione) e con strutture polifase riduce molto il costo.Elaborazione multirate - decimazione e interpolazione →, Progetto di filtri FIR con il metodo delle finestreUn filtro FIR a fase lineare di tipo I (ordine $N$ pari, $h[n]=h[N-n]$) ha risposta $H(e^{j\hat\omega})=e^{-j\hat\omega N/2}\bar H(\hat\omega)$ con ampiezza $\bar H(\hat\omega)=\sum_{n=0}^{N/2}p_n\cos(n\hat\omega)$. Per approssimare un'ampiezza desiderata $D(\hat\omega)$ (passa-basso: $1$ in banda passante, $0$ in banda oscura, con tolleranze $\delta_p,\delta_s$ e frequenze $\hat\omega_p,\hat\omega_s$) ci sono tre metodi. Finestre: si tronca la serie di Fourier di $D$, cioè $h[n]=h_d[n-N/2],w[n]$ con $h_d[n]=\frac{\hat\omega_0}{\pi}\operatorname{sinc}\frac{\hat\omega_0 n}{\pi}$; la rettangolare è ottima in errore quadratico ma dà il fenomeno di Gibbs (9%, $21$ dB), le finestre rastremate (Hann, Hamming, Blackman, Kaiser) abbassano i lobi laterali allargando la transizione ($\hat\omega_s-\hat\omega_p\simeq\alpha,2\pi/L$). Campionamento in frequenza: $h=\mathrm{IDFT}$ dei campioni di $D$, esatto solo sui campioni. Minimax (Parks-McClellan): errore pesato minimo nel caso peggiore, soluzione equiripple con almeno $r+2$ alternanze, ordine più basso a parità di specifiche.Progetto di filtri FIR con il metodo delle finestre →, Interpolazione LTI e teorema del campionamentoIl campionatore R→Z(T) ripete lo spettro con periodo Fc = 1/T (Sc(f) = Σ S(f - kFc), senza fattore 1/T). Il filtro interpolatore Z(T)→R ha y(t) = Σ x(nT) g0(t-nT) con g0 = T g e in frequenza Y = G·X. Se S è nulla fuori da (-B,B) e Fc ≥ 2B, con g0(t) = sinc(Fc t) si ricostruisce esattamente s(t) dai campioni. Altrimenti c'è un errore (in banda per l'aliasing, fuori banda per la parte tagliata), ridotto da un prefiltro anti-aliasing.Interpolazione LTI e teorema del campionamento →, Campionamento e ricostruzioneIl campionamento $s_c(nT)=s(nT)$ trasforma un segnale continuo in uno discreto e, in frequenza, ripete lo spettro con periodo $F_c=1/T$: $S_c(f)=\sum_kS(f-kF_c)$. Se le repliche si sovrappongono si ha aliasing e il segnale non è recuperabile. Un interpolatore $\mathbb Z(T)\to\mathbb R$ con risposta impulsiva $g$ produce $\tilde s(t)=\sum_nT g(t-nT)s(nT)$ e $\tilde S=G,S_c$. Teorema del campionamento: se $s$ ha banda $B$ e $F_c\ge2B$, l'interpolatore ideale $G=\operatorname{rect}(f/F_c)$ ricostruisce esattamente $s(t)=\sum s(nT)\operatorname{sinc}(F_c(t-nT))$. Se le ipotesi non valgono c'è un errore, in banda e fuori banda, riducibile con un prefiltro anti-aliasing.Campionamento e ricostruzione →.

Codice

import numpy as np
from PIL import Image
from scipy import signal, ndimage
import matplotlib.pyplot as plt

# ---------- 1) Interpolazione lineare di una sequenza 1D, L = 3 ----------
rng = np.random.default_rng(1)
samp = np.round(rng.uniform(-1, 1, 8), 2)            # 8 campioni casuali
L = 3
xss = np.zeros(L * len(samp)); xss[::L] = samp        # espansore: L-1 zeri tra due campioni
h = np.r_[np.arange(1, L + 1), np.arange(L - 1, 0, -1)] / L   # triangolo di picco 1, lunghezza 2L-1
y_full = np.convolve(xss, h)                          # convoluzione completa
y_same = np.convolve(xss, h, mode='same')             # stessa lunghezza dell'ingresso
print("campioni      :", samp)
print("h             :", np.round(h, 4))
print("full:  y[n] = x[n-(L-1)]? ", np.allclose(y_full[L - 1::L][:len(samp)], samp))
print("same:  y[n] = x[n]?       ", np.allclose(y_same[::L], samp))
print("y_full[:12]   :", np.round(y_full[:12], 3))

# ---------- 2) Zoom di una patch di Lena.tiff ----------
img = np.asarray(Image.open('Lena.tiff').convert('L')).astype(np.float64) / 255   # 512x512 in [0,1]
patch = img[240:290, 240:290]                         # patch 50x50 con un occhio

tri = np.r_[np.arange(1, 6), np.arange(4, 0, -1)] / 5           # triangolo, L = 5
k = np.arange(-10, 11); sinc = np.sinc(k / 5)                   # sinc troncato a 21 coefficienti
sinc_h = sinc * np.hamming(21)                                  # sinc finestrato (Hamming di 21 punti, centrata)
print("guadagno in continua (somma dei coefficienti): triangolo", tri.sum(), " sinc", round(sinc.sum(), 3), " sinc+Hamming", round(sinc_h.sum(), 3))

def filtra(z, h1d):                                   # filtro separabile: prima le righe, poi le colonne
    z = signal.convolve2d(z, h1d[None, :], mode='same')
    return signal.convolve2d(z, h1d[:, None], mode='same')

def zoom(low):                                        # i quattro metodi, fattore 5
    zer = np.zeros((5 * low.shape[0], 5 * low.shape[1])); zer[::5, ::5] = low    # inserzione degli zeri
    return {"ripetizione": np.repeat(np.repeat(low, 5, axis=0), 5, axis=1),
            "lineare": filtra(zer, tri),
            "sinc rettangolare": filtra(zer, sinc),
            "sinc + Hamming": filtra(zer, sinc_h)}

for nome, z in zoom(patch).items():                   # patch 50x50 -> 250x250
    print(f"{nome:18s} media interna = {z[15:-15, 15:-15].mean():.4f}   errore sui nodi = {np.abs(z[::5, ::5] - patch).max():.1e}")

# ---------- misura dell'errore: serve una "verita'" ad alta risoluzione ----------
def psnr(a, b, m=15):                                 # PSNR sull'interno (si scartano i bordi)
    a, b = a[m:-m, m:-m], b[m:-m, m:-m]
    return 10 * np.log10(1 / np.mean((a - b) ** 2))

crop = img[140:390, 140:390]                          # 250x250 attorno all'occhio (contiene la patch)
hi = ndimage.gaussian_filter(crop, 1.5)               # riferimento limitato in banda
low = hi[::5, ::5]                                    # 50x50 da ingrandire 5 volte
for nome, z in zoom(low).items():
    print(f"{nome:18s} PSNR = {psnr(z, hi):5.2f} dB")
low_grezza = crop[::5, ::5]                           # stessa prova senza filtro prima della decimazione
for nome, z in zoom(low_grezza).items():
    print(f"{nome:18s} PSNR (rispetto all'originale) = {psnr(z, crop):5.2f} dB")

# ---------- 3) Risposte in frequenza dei filtri 1D ----------
for nome, f in (("lineare", tri), ("sinc", sinc), ("sinc + Hamming", sinc_h)):
    w, H = signal.freqz(f, 1, worN=4096)
    mag = np.abs(H) / 5                                # guadagno normalizzato a 1
    punto = lambda x: mag[np.argmin(np.abs(w - x * np.pi))]
    print(f"{nome:15s} |H|/5 a 0, 0.1pi, 0.2pi, 0.3pi: {punto(0):.3f} {punto(0.1):.3f} {punto(0.2):.3f} {punto(0.3):.3f};"
          f"  massimo oltre 0.4pi: {mag[w >= 0.4 * np.pi].max():.4f}")
    plt.plot(w / np.pi, mag, label=nome)
plt.xlabel(r"$\hat\omega/\pi$"); plt.legend(); plt.grid(); plt.savefig("lab4_risposte.png")

# ---------- 4) confronto con OpenCV ----------
import cv2
cv_lin = cv2.resize(low, (250, 250), interpolation=cv2.INTER_LINEAR)
cv_cub = cv2.resize(low, (250, 250), interpolation=cv2.INTER_CUBIC)
print("cv2 lineare, PSNR", psnr(cv_lin, hi), "-> con spostamento di 2 pixel", psnr(cv_lin[2:, 2:], hi[:-2, :-2]))

Risultati e spiegazione

Riscaldamento. Con $L=3$ il triangolo è $h=\{\frac13,\frac23,1,\frac23,\frac13\}$ (picco $1$, lunghezza $2L-1=5$, somma $L=3$). La convoluzione completa np.convolve(xss, h) restituisce l'ingresso ritardato di $L-1=2$ campioni: $y[nL+(L-1)]=x[n]$ (il filtro è non causale se centrato, ma convolve lo tratta come causale con il picco al centro dei $5$ coefficienti, cioè in posizione $2$); l'output nei punti $n=2,5,8,\dots$ è $0{,}02;\,0{,}9;\,-0{,}71;\dots$ come i campioni originali. Con mode='same' l'uscita viene riallineata e $y[nL]=x[n]$ (verificato). Fra due campioni l'uscita segue i segmenti: dopo $0{,}02$ e $0{,}9$ si hanno $0{,}313$ e $0{,}607$. La proprietà $y[nL]=x[n]$ vale perché $h$ si annulla nei multipli non nulli di $L$ (interpolazione corretta).

Zoom della patch di Lena. Con la patch della consegna ($50\times50$, livello medio $0{,}513$) i quattro metodi danno immagini $250\times250$ e queste proprietà:

(La media è calcolata scartando $15$ pixel di bordo, dove la convoluzione con zeri falsa il risultato; è confrontabile con quella della ripetizione, $0{,}4868$.)

• Il sinc troncato sembra il metodo ideale ma ha guadagno in continua sbagliato: la somma dei $21$ coefficienti è $4{,}531$ invece di $5$ (la serie $\sum\operatorname{sinc}(k/5)=5$ converge lentamente, e il troncamento taglia una coda che si somma con segni alternati). In 2D il guadagno è al quadrato: $(4{,}531/5)^2=0{,}821$, e infatti la media dell'immagine ingrandita è $0{,}3988/0{,}4868=0{,}819$ di quella della ripetizione: l'immagine è più scura del $18\%$ circa. In più il troncamento dà ondulazioni di Gibbs attorno ai bordi (nei filtri, il sovraelongo di $1{,}12$ nella banda passante, vedi sotto). La finestra di Hamming porta la somma a $5{,}021$ (guadagno 2D $1{,}008$) e abbassa i lobi laterali.
• La ripetizione non aggiunge informazione e dà l'aspetto a blocchi (i pixel replicati sono un'interpolazione "di ordine zero").
• In tutti e quattro i casi i pixel originali sono conservati nei nodi (errore $\lesssim10^{-16}$): tutti i filtri hanno zeri nei multipli di $5$ (interpolazione corretta).

Quanto è buono lo zoom? Sulla patch non c'è una "verità" con cui confrontarsi: l'unica immagine disponibile è quella piccola. Per misurare l'errore si prende un'area più grande di Lena ($250\times250$ pixel, img[140:390, 140:390], centrata sull'occhio), si decima di $5$ per ottenere la patch $50\times50$ e si confronta l'ingrandimento con l'area di partenza, con il PSNR (rapporto tra il quadrato del massimo, $1$, e l'errore quadratico medio, in dB; interno senza $15$ pixel di bordo):

• Area limitata in banda (filtro gaussiano con $\sigma=1{,}5$ pixel prima della decimazione, così non c'è aliasing e il segnale è compatibile con l'interpolazione a banda $\frac\pi5$): il migliore è sinc + Hamming, $1{,}3$ dB sopra il lineare; il sinc troncato è peggiore anche della ripetizione, perché perde il $18\%$ di luminosità (MSE, errore quadratico medio, $0{,}0095$ contro $0{,}0058$ della ripetizione). Su altre due aree di Lena (angolo in alto a sinistra e basso a destra) l'ordine è lo stesso: il lineare dà $37{,}0$ e $34{,}1$ dB, sinc + Hamming $37{,}7$ e $35{,}5$ dB, cioè da $0{,}7$ a $1{,}4$ dB in più.
• Area originale, decimata senza prefiltro: i dettagli (ciglia, bordo del cappello) hanno componenti oltre $\frac\pi5$ che la decimazione ha trasformato in aliasing (Elaborazione multirate - decimazione e interpolazioneL'elaborazione multirate cambia la frequenza di campionamento di un segnale. Interpolazione per $L$ (da $F_s$ a $LF_s$): l'espansore inserisce $L-1$ zeri fra i campioni, $W(\hat\omega)=X(L\hat\omega)$, lo spettro si ripete $L$ volte nell'asse più ampio; un passabasso con taglio $\pi/L$ e guadagno $L$ elimina le immagini, e con il sinc ideale $h[n]=\operatorname{sinc}(n/L)$ i campioni originali sono conservati ($y[nL]=x[n]$); in pratica FIR a finestre o interpolatore lineare (triangolo). Decimazione per $M$ (da $F_s$ a $F_s/M$): $y[n]=x[nM]$, $Y(\hat\omega)=\frac1M\sum_{k=0}^{M-1}X\big(\frac{\hat\omega-2\pi k}M\big)$, repliche sovrapposte (aliasing) se la banda supera $\pi/M$: si fa precedere da un passabasso con taglio $\pi/M$. Conversione razionale $F_s'=\frac LMF_s$ ($L,M$ coprimi): espansore $L$, un solo passabasso di taglio $\min(\frac\pi L,\frac\pi M)$ e guadagno $L$, decimatore $M$; $y[n]=\sum_kx[k],g[nM-kL]$ non è una convoluzione e il sistema è periodicamente tempo-invariante. Realizzare il tutto in più stadi (fattori piccoli prima per l'interpolazione, grandi prima per la decimazione) e con strutture polifase riduce molto il costo.Elaborazione multirate - decimazione e interpolazione →), e nessun interpolatore può recuperarle: il PSNR scende a circa $24$ dB per i due filtri buoni, che diventano praticamente equivalenti ($24{,}07$ e $24{,}02$ dB). L'ingrandimento non può essere meglio di quanto l'immagine piccola permetta.
• Il triangolo dà un risultato buono con soli $9$ coefficienti non nulli contro i $21$ del sinc, per $1{,}3$ dB in meno.

Risposte in frequenza (modulo normalizzato a $1$, cioè diviso per $L=5$; non dipendono dall'immagine):

Il filtro ideale avrebbe $1$ fino a $\frac\pi5=0{,}2\pi$ e $0$ oltre. Il lineare comincia a calare già da $0{,}1\pi$ (attenua le frequenze utili: immagine meno nitida) e non taglia bene le immagini; il sinc troncato ha guadagno $0{,}906$ in continua (da cui il $18\%$ di luminosità persa in 2D) e un sovraelongo in banda passante ($1{,}12$, onda di Gibbs) con lobi laterali moderati; sinc + Hamming ha banda passante piatta (circa $1$) e banda oscura molto bassa.

Confronto con cv2.resize. In OpenCV l'ingrandimento predefinito INTER_LINEAR è un'interpolazione bilineare: equivale al triangolo di questo laboratorio, ma con un allineamento diverso: OpenCV centra i pixel a mezzo campione, cioè il pixel $i$ della patch finisce nella posizione $5i+2$ e non in $5i$ (l'errore sui nodi è $0{,}22$ invece di $0$). Sull'area di prova il PSNR misurato senza correggere lo spostamento è $24{,}00$ dB; spostando di $2$ pixel il confronto (cv_lin[2:, 2:] contro hi[:-2, :-2]) sale a $31{,}66$ dB, praticamente il valore del triangolo ($31{,}59$ dB): sono lo stesso filtro. L'opzione INTER_CUBIC (interpolazione cubica, non è un sinc) allineata arriva a $33{,}62$ dB, più di sinc + Hamming ($32{,}93$ dB): un filtro di soli $4$ punti per direzione ben progettato batte un sinc a $21$ coefficienti.