Algoritmi per la compressione di immagini a istogramma sparso nella catena jpeg2000

CAPITOLO 1. WAVELET E MULTIRISOLUZIONE 13
 10  8  6  4  2 0 2 4 6 8 10
 1
 0.8
 0.6
 0.4
 0.2
0
0.2
0.4
0.6
0.8
1
Figura 1.1: Rappresentazione di Fourier di un'onda quadra
Figura 1.2: Un esempio di Waveletdominio del tempo a quello della frequenza comporta la perdita totale d'informazionespaziale ora, sostituendo il concetto di frequenza con quello di scala (Figura 1.3) invece
Figura 1.3: Il segnale visto dalla trasformata Waveletdi spostarsi da una frequenza ad un'altra, ci si sposta da una risoluzione piu ne a unapiu grossolana. Si colgono cos in un colpo solo sia i dettagli sia le informazioni chepotrebbero essere denite \di massima". Il fatto che nel graco di Figura 1.3 le aree
CAPITOLO 1. WAVELET E MULTIRISOLUZIONE 14tracciate abbiano la stessa supercie e l'esemplicazione del tentativo della trasformataWavelet di raggirare l'imposizione del principio di Heisenberg [26]:Non e possibile determinare esattamente e contemporaneamente posizione equantita di moto di un corpuscolo.Riportato in questo campo signica che non e possibile avere pari accuratezza dal puntodi vista del tempo e della frequenza nell'analisi di un segnale.Per dare un'idea di come lavora la trasformata si prenda il segnale x(t)=sin(t)einun punto a caso si inserisca una discontinuita impercettibile (per esempio a t = ). Latrasformata di Fourier di questo segnale mostra solo un picco sulla frequenza interessata eanch'essa, come l'occhio umano, non riesce a evidenziare la piccola discontinuita (Figura1.4).
Figura 1.4: Una sinusoide imperfetta e la sua trasformata di FourierUsando invece la trasformata Wavelet e possibile trattare separatamente i dettaglie le approssimazioni del segnale arrivando ad un'immagine di questo tipo (Figura 1.5).L'immagine delle approssimazioni (sulla sinistra) corrisponde alla visione denita \dimassima" e infatti anche qui non vi e traccia dell'imperfezione (questa e una caratteristi-ca utile quando per esempio si debbano analizzare dei trend), ma osservando i coeÆcientidei dettagli (sulla destra) si nota che i valori in modulo maggiore sono posti esattamentenella posizione dove la discontinuita era stata in origine inserita, mantenendo quindianche l'informazione temporale.
CAPITOLO 1. WAVELET E MULTIRISOLUZIONE 15
0 100 200 300 400
 1.5
 1
 0.5
0
0.5
1
1.5
0 100 200 300 400
0
1
2
3
4
5
6
x 10
 3
Figura 1.5: Approssimazioni e dettagli della stessa sinusoideNon occorre molto per accorgersi della potenza di questo strumento: la lista chesegue abbraccia solo una parte delle possibili applicazioni: Nel campo monodimensionale (segnali da R a R){ eliminazione del fruscio dai le musicali ([3], [16]){ studio di onde sismiche [37]{ studio dei trend nanziari [20] Nel campo bidimensionale (segnali da R2 a R){ svariate applicazioni in campo medico ([7], [11]) e{ tutto cioche concerne l'elaborazione d'immagine quindi rilevamento contornie caratteristiche, compressione. . . [19]1.2 La Trasformata Wavelet ContinuaLe prossime tre sezioni fanno da sfondo teorico alla comprensione della trasformataWavelet: dopo averne dato una denizione e illustrato le proprieta di base, si passa avedere la trasformata discreta eettivamente usata nell'elaborazione di immagini.
CAPITOLO 1. WAVELET E MULTIRISOLUZIONE 161.2.1 DenizioneUn approccio valido non solo nel campo delle Wavelet e quello di modellare il segnalecome una combinazione lineare di funzioni: per rappresentare il vettore v si consideriuna base b per lo spazio vettoriale V, quindi8v 2 V v = nXk=1 kbk (1.1)dove k e il k-esimo coeÆciente e bk e il k-esimo vettore della base. Passando allefunzioni si puo usare lo stesso approccio:f(t)= nXk=1 kk(t) (1.2)dove in questo caso k(t) e la k-esima funzione di base; se poi f e g soddisfano lacondizione Z +11 jf j2 < +1 (1.3)ovvero f 2 L2(R), si puo denire il prodotto scalare tra due funzioni (
 denota ilcomplesso coniugato): <f;g>= Z ba f(t)g(t)dt (1.4)A questo punto la denizione di CWT1 si vede come prodotto scalare tra la MotherWavelet  e il segnale x(t) usando la (1.4):CWT x (;s)=	 x (;s)=Z ba x(t) ;s(t)dt (1.5)dove  ;s = 1ps   t s  (1.6)La 1.6 deve godere di una proprieta importante, chiamata ipotesi di regolarita:Z ba  dt =0 (1.7)cioe la Wavelet generatrice deve avere valor medio nullo nel dominio del tempo, quindiun andamento oscillatorio (dall'inglese wave, onda).1Continous Wavelet Transform, trasformata Wavelet continua.
CAPITOLO 1. WAVELET E MULTIRISOLUZIONE 17Il momento della Wavelet Espandendo la (1.5) in serie di Taylor no al n-esimoordine centrata in t =0e =0,otteniamo:T (s; 0) = 1ps " nXk=0 x(k)(0) Z tkk!   ts dt+O(n+1)# (1.8)Se a questo punto deniamo il momento della Wavelet Mk comeMk = Z tk (t)dt (1.9)possiamo riscrivere la (1.8) come1ps " nXk=0 x(k)(0)k! Mk(t)sk+1 + O(sn+2)# (1.10)Si sa che il primo termine della sommatoria di (1.10) deve essere nullo per la condizionedi regolarita (1.7). Se si fa in modo anche gli altri termini si annullano la funzionedecadra a zero come s(n+2) e la trasformata Wavelet sara di ordine n.1.2.2 Ortogonalita e ortonormalitaDue vettori v e w si dicono ortogonali se il loro prodotto scalare enullo, in altri termini:< v;w >= nXk=1 vk 
 wk =0 (1.11)e un insieme di vettori fv1;v2;:::;vng si dice ortonormale se vale anche:8v jvj2 =1 (1.12)Similmente due funzioni f e g sono ortogonali se si annulla la (1.4) e un insieme difunzioni f1;2;:::;ng e ortonormale se inoltreZ ba jk(t)j2dt =1 (1.13)La proprieta di ortonormalitae utile nella computazione delle Wavelet perche permettedi calcolare il coeÆciente k in maniera veloce:k =<f;k >= Z ba fkdt (1.14)
CAPITOLO 1. WAVELET E MULTIRISOLUZIONE 18quindi se vale la (1.2), allora f = nXk=1 <f;k >k (1.15)Non sempre queste proprietasonodisponibili, quindi si deve ricorrere ai concetti menoforti di biortogonalita (quest'ultima implica che si usino due basi ortogonali tra di loroma non formanti una base ortonormale) e frame: condizione meno forte della 1.15, checomunque fornisce un limite inferiore e superiore a f :A k f k2 Z ba (f ;s)2  B k f k2 A;B 2 R (1.16)La teoria necessaria ad un uso pratico delle Wavelet e stata sviluppata da Mal-lat [18], che ha introdotto l'analisi in multirisoluzione. Supposto che f soddis la 1.3,viene denito un operatore A2j che realizza l'approssimazione di f(x) alla risoluzione2j, denita come A2j 
 f(x), che puo essere vista come la proiezione ortogonale di f(x)sullo spazio V2j  L2(R). Lo spazio V2j puo essere denito comef(
) 2 V2j () f(
)e costante in  k2j ; k+12j  (1.17)Una base creata con funzioni cos denite e soddisfacente solo a livello teorico in quantoqueste non godono di alcun tipo di regolarita.Per caratterizzare numericamente l'operatore A2j bisogna trovare una base ortonor-male per lo spazio V2j . Il teorema dimostrato in [17] assicura che9!(x):8j 2 Z2j(x)=2j(2jx)=) 1p2j 2j (x n2j)eunabase ortonormale per V2j(1.18)La funzione (
) viene chiamata scaling function e funge da ltro passa basso per ilsegnale f(x).1.2.3 Dalle Wavelet continue a quelle discretePrimo problema: discretizzazione delle scalature e delle traslazioniLa prima diÆcolta computazionale delle Wavelet e quella di avere un parametro s discalatura e  di traslazione temporale continui: per ovviare a questo problema si riscrive
CAPITOLO 1. WAVELET E MULTIRISOLUZIONE 19
la (1.6) in maniera discreta: j;k(t)= 1qs j0   t k0s j0s j0 ! (1.19)dove j; k 2 N e s0 > 1e il passo ssato per lo scaling. Il passo di traslazione 0 e funzionedella scalatura. I passi comuni sono s0 =2e0 = 1. Dopo questa discretizzazione, se valela proprieta 1.16, il segnale originario puo essere ricostruito a partire dai coeÆcienti delleWavelet. In questo caso il termine  j;k viene chiamato frame con limiti A e B.SeA e Bsono uguali, il frame viene detto tight, altrimenti bisogna ricorrere per la ricostruzionedel segnale al dual frame. Dal momento che le Wavelet sono discrete, possono essererese ortonormali facendo in modo che rispettino la (1.13); in questo modo il segnale puoessere ricostruito applicando il concetto espresso in (1.2), che porta a:f(t)=Xj;k 
(j; k) j;k(t) (1.20)Secondo problema: limitazione della bandaAnche discretizzando la scalatura e la traslazione si ha comunque una componente in-nita nelle frequenze: il problema puo essere raggirato osservando che la Wavelet hauno spettro di tipo passa-banda, quindi, per ricoprire tutte le frequenze, si puo usare laproprieta di scalatura che viene dalla trasformata di Fourier:Fff(at)g = 1jajF !a (1.21)Comprimendo nel tempo la Mother Wavelet si opera un allargamento dello spettro dellatrasformata. Rimane il problema di come coprire l'ultimo tratto di frequenza che si portano alla componente continua, visto che dimezzando sempre la banda non si potra maigiungere a frequenza nulla. In questo caso si ricorre alla scaling function (introdotta in1.18), un ltro passa basso che ricopre proprio l'area di interesse.Terzo problema: come calcolare la trasformataLa tecnica usata per risolvere il terzo ed ultimo problema echiamata codica di sottoban-da (introdotta a livello teorico nella Sezione precedente): si riprende in toto il concetto
CAPITOLO 1. WAVELET E MULTIRISOLUZIONE 20
Figura 1.6: La compressione temporale sulla Mother Wavelet produce una dilatazionein frequenza.
Figura 1.7: La funzione di scaling funge da ltro passa basso per terminare la coperturain frequenza delle altre Waveletdei ltri passa banda di Figura 1.6, e si itera il processo no a quando si arrivaadunalarghezza di banda soddisfacente al livello di dettaglio che si vuole ottenere (per ognidecomposizione, il ltro passa basso corrispondente fornisce invece un'approssimazionedella funzione d'ingresso).1.2.4 Implementazione della trasformata Wavelet discretaApplicazione ricorsiva dei ltriSupponendo che h[n] sia la risposta all'impulso del ltro passa basso mentre g[n] quelladel passa banda, il ltraggio risulta uguale ad una convoluzione del segnale con questeultime: x[n] ?h[n]= 1Xk=1x[k] h[n k] (1.22)
CAPITOLO 1. WAVELET E MULTIRISOLUZIONE 21Questa operazione, secondo il teorema di Nyquist, permette di dimezzare il numero dicampioni del segnale, visto che l'uscita del passa banda e del passa basso hanno ridottola frequenza di un fattore 2. Questa eliminazione prende il nome di sottocampionamentoe puo essere implementata insieme alla convoluzione in questo modo:y[n]= 1Xk=1h[k] x[2n k] (1.23)Basandosi sul concetto espresso nella formula (1.22) e usando la (1.23), la trasformataWavelet discreta (la decomposizione di primo livello) viene cos implementata:yH [k] = Xn x[n] g[2k  n] (1.24)yL[k] = Xn x[n] h[2k  n] (1.25)Supponendo che il segnale sia campionato ad una frequenza di 2 rad/s, quindi lospettro di frequenza vada da 0 a  rad/s, e abbia 1024 campioni, dopo il primo ltraggiole uscite dai due ltri avranno 512 campioni a testa, ma il primo segnale avra uno spettronell'arco [0;=2] rad/s e il secondo in [=2;] rad/s. L'operazione viene ripetuta e i 4segnali risultanti copriranno rispettivamente le frequenze [0;=4], [=4;=2], [=2; 3=4]e[3=4;] rad/s, ma avendo 256 campioni ciascuno. Ad ogni passaggio viene raddoppiatala frequenza e dimezzata la risoluzione: in questo caso si ripete il procedimento per altre7 volte no ad arrivare a segnali con 2 campioni ciascuno. Il segnale viene ricostruitoconcatenando i coeÆcienti dalla coda alla testa della catena di ltraggi, facendo cosinmodo che il segnale e la DWT abbiano lo stesso numero di coeÆcienti. Il ltro passabasso e il ltro passa banda non sono senza alcun legame, ma la loro relazione e:g[L 1 n]=(1)n h[n] (1.26)dove L e il numero di punti del ltro. La formula per ricostruire il segnale e moltosemplice, dal momento che i ltri formano una base ortonormale:x[n]= 1Xk=1 (yH [k] g[2k  n])+(yL[k] h[2k  n]) (1.27)
Capitolo 2Stato dell'arte
In questo capitolo viene presentato lo stato dell'arte nell'ambito della compressione delleimmagini sse precedente al JPEG2000 ([1], [27]) sia nella compressione lossless, conJPEG-LS [33], sia nella compressione lossy (a perdita d'informazione) con l'algoritmoEZW [25]. Nella prima sezione viene presentato anche il rapporto con la sparsita delleimmagini, che verra ripreso poi nel Capitolo 5; la seconda sezione si concentra invecesulla bonta delle Wavelet nella compressione delle immagini.2.1 Istogrammi sparsi e compressione lossless2.1.1 IntroduzioneDa diversi campi scientici quali l'acquisizione di immagini mediche o spaziali la richiestadi compressione senza perdita di dati e sempre piu pressante. In questo tipo di immaginie importante anche l'informazione contenuta nel singolo pixel, quindi non e possibileperdere nessun tipo di informazione con la compressione; nel 1999 l'organizzazione ISOe arrivata alla denizione dello standard JPEG-LS1.Oltre al JPEG-LS, altri due algoritmi si sono aacciati sulla scena informatica:CALIC [36] e JPEG2000 : nel primo caso i risultati raggiunti erano molto miglioririspetto a quelli del JPEG-LS ma la complessita computazionale richiesta lo rendeva1LosslesS 22
CAPITOLO 2. STATO DELL'ARTE 23non competitivo, nel secondo caso l'uso della trasformata Wavelet discreta anziche lapredizione usata nel JPEG-LS migliora di molto le prestazioni senza un eccessivo sprecodi risorse.In tutti e tre le situazioni viste sopra, l'uso di immagini con istogramma sparsopeggiora di molto i rapporti di compressione. In questa sezione si esamina questo fattorapportato ai tre diversi formati e si presenta un algoritmo di preprocessing on-line dainserire nella catena JPEG-LS per cercare di ovviare a questo problema.2.1.2 Come lavora il JPEG-LSIl JPEG-LS puo essere analizzato dividendolo in tre blocchi principali: predizione, mo-dellizzazione e codica, che verranno illustrate nelle sezioni successive.PredizioneChiamato x il campione corrente dell'immagine, in Figura 2.1 sono mostrati i quattrocampioni che vengono usati per la predizione. La parte ssa del predittore assegna al
Figura 2.1: Le vicinanze del campione xcampione successivo il valore x^t+1 in base a questo schema:x^t+1  8>><>>: min(a; b) se c  max(a; b)max(a; b) se c  min(a; b)a + b c; altrimenti (2.1)Questo predittore puo anche essere visto come un rilevatore di contorni mediano chesceglie tra i tre predittori ssi a, b, a + b  c. Il campione d non compare nella (2.1)perche viene usato solo nella parte adattativa del predittore.
CAPITOLO 2. STATO DELL'ARTE 24Modellizzazione del contestoIn questa parte si utilizza il campione d non preso prima in considerazione: vengonocalcolate le dierenze g1 = d bg2 = b cg3 = c a (2.2)Per ridurre il numero di informazioni derivanti dalla 2.2, le coppie di dierenze gi;jvengono processate con un quantizzatore a step equiprobabili; a questo punto si calcolala probabilita condizionale tra la dierenza di x^t+1 e xt+1 e il modello delle triplette didierenze. Il modello di probabilitache viene usato prende il nome di TGSD2:P(;s)()= (1 )1s + s j+sj (2.3)dove  e la dierenza tra il campione predetto e quello eettivo mentre , compreso tra0 e 1, controlla il decadimento esponenziale della curva e s e un parametro usato perbilanciare l'oset introdotto nella fase di predizione ssa. L'andamento della curva eillustrato in Figura 2.2.
 3  2  1 0 1 2 3
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
epsilon
P
s 
Figura 2.2: La distribuzione geometrica a due code2Two Sided Geometric Distribution