[ITA] Digital image compression tecniques

Compressione di immagini digitali. Algoritmi a confronto. Andrea Favilli 6 Marzo 2014 Ridurre le dimensioni delle immagini digitali penalizzando il meno possibile la qualit` a visiva: questa l’ambizione di ci`o che andremo a trattare. Vedremo i principi che stanno dietro uno degli algoritmi pi` u diffusi, il JPEG, acronimo di Joint Photographic Experts Group, con relativa implementazione e confronto con le possibili alternative. Capiremo quali sono gli effettivi vantaggi che portano tale standard ad una cos`ı ampia diffusione. Vista l’esigenza di una computazione efficiente e rapida andremo ad utilizzare il pacchetto commerciale Matlab, in luogo dell’alternativa gratuita GNU Octave: il fatto che Matlab abbia a disposizione molti degli algoritmi matematici indispensabili nella nostra trattazione ci ha fatto scartare l’uso di un linguaggio a pi` u basso livello di astrazione, penalizzando l’efficacia finale, ma rimanendo nel puro spirito sperimentale. Useremo shell scripting da terminale Linux per coordinare le operazioni delle singole funzioni e piccole parti di C per poter gestire i file in maniera agevole. Gli algoritmi in oggetto sono essenzialmente tre: il filo conduttore della nostra esperienza sta nel fatto che non opereremo “una tantum” sulla singola immagine, ma al contrario suddivideremo questa in sottounit`a (di dimensione variabile) e su tali oggetti andremo ad effettuare compressione. Non perdendoci in ulteriori spiegazioni anticipate, iniziamo immediatamente con la trattazione.

Alcuni prerequisiti: codifiche di immagini. L’immagine digitale `e concepita come un insieme di unit`a minime (dei piccoli quadratini) chiamate pixel, ciascuno di essi pu`o assumere una certa luminosit`a e uno specifico tono di colore. La configurazione dei pixel `e regolata da una o pi` u matrici scritte secondo opportune regole di codifica: il fatto che un’immagine sia esprimibile in forma matriciale fa da subito riflettere su come molte delle conoscenze messe a disposizione dall’algebra lineare siano utili al trattamento della grafica. Vediamo rapidamente come sono fatti alcuni file: esaminiamo quei formati ad alto dispendio di spazio, costruiti a partire dalle risposte del sensore luminoso delle fotocamere. In poche parole, questo tipo di foto hanno una qualit`a pressoch`e assimilabile a quella dell’analogico, ma non sono per l’appunto gestibili. Le immagini in scala di grigio sono codificate nel formato .pgm (portable gray map): `e possibile costruire un file di questo tipo prendendo un editor di testo e inserendo determinate informazioni secondo uno specifico ordine. 1

• Un codice identificativo del tipo di file, nel caso del .pgm scriveremo P2. • Eventuali commenti preceduti dal carattere #. • Le dimensioni della foto, rispettivamente larghezza e altezza. • La massima luminosit` a che il singolo pixel pu`o assumere, lo standard `e 255. • I fattori di luminosit` a dei singoli pixel (valori numerici da 0, il nero, a 255, il bianco). I pixel vengono esplorati per riga da sinistra verso destra. Riportiamo di seguito un possibile esempio: P2 #Immagine di prova. 1024 768 255 198 112 112 113 115 ... Da segnalare il fatto che i fattori di luminosit`a sono talvolta espressi in forma binaria e dunque non leggibili con un comune editor. Fortunatamente, alcuni software per la grafica (tra cui includiamo il blasonato Gimp) permettono la scelta della formattazione (binaria o ASCII) al momento del salvataggio dell’immagine. Un estensione di ci` o che si `e appena mostrato `e il file .ppm (portable pix map), volto a codificare immagini a colori. Tale formato `e caratterizzao da: • Il codice identificativo P3. • Eventuali commenti preceduti da #. • Le dimensioni dell’immagine, rispettivamente larghezza e altezza. • Il massimo indice di intensit`a colore, di default 255. • Le triple recanti la configurazione di ciascun pixel secondo RGB (vedi paragrafo successivo), i pixel sono esplorati anche in questo caso per riga da sinistra verso destra. I due esempi di codifica che abbiamo visto sono alcuni dei possibili, ma sono quelli che sono utili al nostro scopo. Esistono numerose estensioni di file immagine non compressi (.pnm, .pbm), oltre ad alcune estensioni proprietarie dei produttori di fotocamere (.cr2 per Canon, .arw per Sony) non facilmente leggibili poich`e codificate secondo regole RAW (binarie) e secondo criteri per l’appunto non standard. Alcune informazioni sui formati canonici possono essere reperite sul sito http: //netpbm.sourceforge.net/doc/#formats.

2

RGB e YCbCr Come posso esprimere la tonalit`a di colore che il singolo pixel pu`o assumere? I modelli (detti ”spazi colore”) sono molteplici, scelti soprattutto in funzione dell’applicazione in questione. In questo paragrafo ne presentiamo due, uno necessario alla codifica del .pgm, mentre l’altro utilizzato nell’ambito del JPEG. • Decompongo il colore in un sistema di toni fondamentali, ossia per ciascuna unit` a ho una tripla di valori corrispondenti alle quantit`a di rosso di verde e di blu. Tale sistema `e per l’appunto detto RGB, da red green blue. In termini complessivi posso pensare ad una foto come alla sovrapposizione di tre immagini monocromatiche: l’una in scala di rossi, la seconda in scala di verdi e la terza in scala di blu. Questo formato `e intuitivamente il pi` u naturale, ma allo stesso tempo non `e molto conveniente dal punto di vista dell’intepretazione da parte di dispositivi elettronici. • Posso invece pensare il pixel rappresentato da una componenente luminosa (detta ”lumia”) con due componenenti cromatiche (dette appunto ”di crominanza”), abbiamo ancora una volta una tripla di valori: questo spazio colore `e detto YCbCr. Una foto `e vista come la combinazione della sua corrispettiva in scala di grigi (Y) e di due foto policromatiche opache (Cb e Cr): le tonalita scure (es. marrone) hanno alti livelli sul canale Cr e bassi su Cb, alcune tonalit` a (come il verde) hanno livelli bassi su entrambi i canali, il Cb invece predilige le tonalit` a accese (es. azzurro cielo). Un vantaggio evidente di tale sistema, oltre ad essere compatibile con il processo di visualizzazione dei monitor, consiste nella separazione della luce dal colore ed `e proprio ci`o che lo rende ottimale per il JPEG. Veniamo al punto fondamentale: come eseguo conversioni dallo spazio colore RGB a quello YCbCr? La trasformazione che associa la tripletta dei colori fondamentali a quella dello spazio colore in lumiocrominanza `e in realt`a un’applicazione lineare. Se R0 , G0 e B 0 sono le componenti RGB shiftate di -128, la tripla Y , Cb, Cr `e data da:      0 Y 0, 299 0, 587 0, 114 R Cb = −0, 1687 −0, 3313 0, 5  G0  0, 5 −0, 4187 −0, 0813 B 0 Cr La conversione dell’YCbCr nella forma RGB segue immediatamente dal calcolo dell’applicazione inversa:  0    1 0 1, 402 R Y G0  = 1 −0, 344 −0, 714 Cb B0 1 1, 772 0 Cr Vediamo un esempio di scomposizione nei due spazi colore che abbiamo esaminato. Il codice utilizzato per ricavare le componenti RGB e YCbCr `e reperibile all’indirizzo: http://poisson.phc.unipi.it/~afavilli/Software/sep_ componenti/.

3

Figura 1: Immagine originale.

Figura 2: Schema RGB.

4

Figura 3: Schema YCbCr.

Alcuni prerequisiti: elementi di shell scripting. Il terminale di Linux `e uno strumento molto potente, in quanto offre sia la possbilit` a di immettere comandi in maniera interattiva che creare particolari file di testo eseguibili (in formato .sh) con i quali eseguiremo in maniera del tutto automatica una sequenza di comandi prestabiliti. Quest’ultima consuetudine `e chiamata shell scripting: tutto ci`o `e assai utile al nostro scopo, creeremo infatti uno script per terminale in grado di far eseguire funzioni Matlab e piccoli programmi in C che nel complesso effettueranno compressione e decompressione di immagini. Osserviamo subito un esempio pratico: creiamo un file run.sh con il suddetto contenuto. #!/bin/sh mkdir prova cd prova touch test Leggiamo riga per riga: la stringa #!/bin/sh indica essenzialmente di eseguire tutto ci` o che segue da shell, ossia creare nella posizione corrente una cartella denominata prova e di spostarsi all’interno di questa per creare il file test. A questo punto apriamo il terminale, impostiamo l’eseguibilit`a del file .sh con chmod -x run.sh ed eseguiamo lo script con sh ./run.sh. Verranno eseguite precisamente le istruzioni proposte. Quella seguita finora `e una procedura standard per la creazione di script: quan5

do nel presente testo lavoreremo con file .sh si assume ogni volta di dover eseguire il chmod dopo il salvataggio. Illustriamo in maniera molto sintetica alcune opzioni che solitamente risultano ridondanti in esecuzione interattiva, ma sono fondamentali per lo scripting. >>> echo "Come ti chiami?" Come ti chiami? >>> read NOME Andrea >>> echo "Ciao, \${NOME}!" Ciao, Andrea! L’output su schermo si ottiene facilmente con la funzione echo, mentre l’input con read seguito dal nome della variabile. Posso richiamare la variabile in qualsiasi punto dello script racchiudendola tra parentesi graffe, la cui aperta `e preceduta dal simbolo $. Ultimo concetto che vediamo `e il passaggio di istruzioni ad una riga di comando esterna, come pu` o essere quella di Matlab o di Gnuplot: basta utilizzare una funzione echo con contentimento delle istruzioni da passare seguita dal simbolo | e dal comando che evoca l’applicazione in oggetto. Un esempio: echo "eig(randn(100))" | matlab Siamo ora pronti ad iniziare con gli algoritmi veri e propri. Ricordo che il codice dell’intera sperimentazione (comprese le variazioni non presenti su questo testo) `e reperibile qui: http://poisson.phc.unipi.it/~afavilli/Software/ Compr_Immagini/.

Primo approccio: una SVD ”a blocchi”. Cerchiamo di ribadire i nostri obiettivi: vogliamo trovare un metodo per poter comprimere una data immagine, cercando il pi` u possibile di mantenere la qualit` a, sapendo che essenzialmente ho in partenza tre matrici. Un’idea ci giunge da alcuni risultati dell’algebra lineare: considero una matrice A ∈ Rm×n , allora esiste sempre una decomposizione (detta ai valori singolari, SVD) nelle matrici U , Σ, V tale che A = U ΣV T dove U ∈ Rm×m , le cui colonne ui sono dette vettori singolari sinistri, V ∈ Rn×n , le cui colonne vi sono dette vettori singolari destri. U e V sono entrambe matrici ortogonali. Σ ∈ Rm×n `e una matrice ”diagonale”, ossia Σij = 0 se i 6= j, Σij = σi se i = j. Gli elementi σi sono detti valori singolari, vale σ1 ≥ . . . ≥ σmin(m,n) . Alla luce di quanto appena detto, possiamo riscrivere la decomposizione in notazione vettoriale: min(m,n) X A= σi ui vTi i=1

6

La propriet` a chiave (dalla quale segue l’algoritmo di compressione) consiste nel fatto che la matrice Ak definita come k X

Ak =

σi ui vTi

i=1

`e quella (di rango k) che minimizza la funzione X 7→ min kA − Xk2 X

ed `e quindi una candidata ideale all’approssimazione di A, essendo consapevoli del fatto che posso calcolare Ak conoscendo solo k valori e k vettori singolari destri e sinistri. Definendo uh,i e vk,i come le coordinate h-esima di ui e k-esima di vi rispettivamente, posso infatti esprimere Ak come (Ak )ij =

k X

σt ui,t vj,t .

t=1

Se A `e una matrice che rappresenta una delle tre componenti di un’immagine (RGB o YCbCr che sia), posso calcolarne la SVD (non parleremo qui dei metodi per farlo, Matlab ci viene in aiuto), conservarne un certo numero k di valori/vettori singolari e codificare questi in un file: questo `e un valido algoritmo di compressione. Per decomprimere mi `e sufficiente ricostruire la matrice ` possibile fissare k  min(m, n), Ak sfruttando la relazione sopra mostrata. E solitamente (per un fissato ε > 0) si sceglie k tale che σσ1i < ε. Per sperimentare `e tuttavia cosa buona lasciare il valore di taglio all’utente in modo da poter studiare bene il calo di qualit` a e di spazio in relazione al k scelto. Non realizzeremo precisamente questo algoritmo, ma vogliamo invece operare con una procedura che sia confrontabile con il JPEG che svilupperemo dopo. Quello che svolgeremo `e una sorta di ibrido: non eseguiremo compressione SVD sull’intera matrice A, ma fisseremo un valore p ∈ N e divideremo la matrice in blocchetti p × p, su ciascuno di questi blocchetti andremo poi ad eseguire compressione come sopra illustrato. Da notare il fatto che le matrici U e V nella decomposizione di ciascun blocco saranno questa volta quadrate e avremo esattamente p valori singolari per blocco. • Un problema si ha qualora su A ∈ Rm×n , m 6= 0 mod p oppure n 6= 0 mod p: affinch`e sia effettivamente possibile la separazione in blocchi equivalenti p × p `e necessario modificare la dimensione di A aggiungendo ˜ Pi` pixel, ottenendo A. u precisamente, se m 6= 0 mod p dovr`o aggiungere p − (m mod p) righe di pixel, se invece n 6= 0 mod p aggiunger`o p − (n mod p) colonne. Per non creare particolari problemi nella compressione andr` o ad assegnare a tali pixels il valore medio tra 0 e 255, ossia 128. Dovr` o in sintesi accertarmi che il mio codice includa una fase iniziale con lo scopo di accertarsi la divisibilit`a in blocchi ed eventualmente creare A˜ come spiegato sopra. Alla fine di questo passaggio ho sicuramente una possibile suddivisione dell’immagine: A˜ = (Bij )ij 7

• Vediamo alcune strategie per risparmiare spazio. Anzitutto `e possibile modificare i vettori singolari in modo da non dover conservareqanche i (ij)

(ij)

(ij)

corrispettivi valori singolari: `e sufficiente definire u ˆ h = uh σh e q (ij) (ij) (ij) m n v ˆ h = vh σh , ∀h = 1, . . . , p, ∀i = 1, . . . , p , ∀j = 1, . . . , p . (Gli (ij) (ij) (ij) u , v ,σ sono i vettori singolari e i valori singolari del sottoblocco Bij ). In tale modo ho ottenuto k X

(ij) (ij) (ij)T

σh uh vh

h=1

=

k X

(ij) (ij)T

u ˆh v ˆh

h=1

Ci` o significa esattamente che posso risparmiare il salvataggio di p valori numerici per blocco, in totale ben ( mn p ) valori! La presenza della radice quadrata nella nostra ridefinizione dei vettori causa per`o una decimalizzazione degli argomenti: per ogni vettore singolare prenderemo soltanto la parte intera delle componenti moltiplicate per un opportuno fattore γ, sceglieremo ad esempio γ = 20, e riusciremo in tal modo a ridurre ancora la quantit` a di memoria occupata. • Una questione che pu` o sorgere naturale giunti a questo punto: la conservazione di k vettori singolari sinistri e k destri rappresenta veramente una riduzione di dati rispetto ad una matrice di p2 elementi? Ciascun vettore singolare ha p componenti, in totale ho 2pk elementi da memorizzare in luogo dei canonici p2 : `e evidente che ci`o `e effettivamente vantaggioso solo se k < p2 . Come vedremo si hanno buoni risultati anche con valori di k effettivamente molto piccoli, quindi questo non rappresenta un effettivo problema. • Lamnmatrice che vado a scrivere su file come immagine compressa `e M ∈ R p ×2k , cos`ı costruita:  (11)  (11) Uk Vk  (21) (21)   Uk  Vk   ..   ..  . .   (12)  (12)  U Vk  k   (22) (22)  Vk  Uk  M = ..   ..  .   .  (13)  (13)  Uk  Vk  (23)  (23)  U Vk  k   . ..   .. .   mn  n (m ( p p) p p) Uk Vk (ij)

(ij)

(ij)

(ij)

(ij)

(ij)

ˆk ). Essenzialmente dove, Uk = (ˆ u1 . . . u ˆ k ) e Vk = (ˆ v1 . . . v vado ad esplorare la matrice A˜ blocco a blocco per colonne, per ciascuna suddivisione vado a calcolare la SVD e accodo i primi k vettori singolari destri e sinistri (opportunamente modificati e rinormalizzati) alla matrice M. 8

Siamo pronti a delineare l’algoritmo per poi giungere alla stesura del codice. 1. Preparazione alla suddivisione in blocchi: eventuale aggiunta di righe e ˜ colonne (da A passo ad A). 2. Esecuzione della compressione: per ogni blocco Bij calcolo [U (ij) , Σ(ij) , U (ij) ] = (ij) (ij) SV D(Bij ), ricavo poi Uk e Vk eseguendo il taglio ai primi k vettori singolari e facendo ”assorbire” i valori singolari agli stessi vettori, eseguo poi la moltiplicazione per γ prendendo solo la parte intera di ogni componente. Aggiorno la matrice M dopo l’iterazione su ogni blocco. 3. Salvataggio della matrice M che rappresenta l’immagine compressa: discuteremo in seguito le modalit`a. (ij)

4. Decompressione: recupero da M i vettori singolari e calcolo (Bk )ab = Pk (ij) (ij) (ij) γ −2 h=1 u ˆ a,h v ˆb,h . Alla fine ottengo A˜k = (Bk )ij . 5. Recupero dimensioni originali: da A˜k ottengo Ak , approssimazione della A di partenza. 6. Riscrittura di A come .ppm. Queste le istruzioni da eseguire con la sola differenza che le matrici dell’immagine a colori saranno tre: dovr`o effettuare la procedura lavorando sulle matrici in parallelo. Prima di vedere il codice `e necessario decidere se implementare il tutto sullo spazio colore RGB o su lumiocrominanza: poich`e le sperimentazioni effettuate non comportano grandi differenze di rendimento e per il fatto che la YCbCr comporta un numero di operazioni totale maggiore opteremo per la RGB. (Il codice analogo in YCbCr `e reperibile all’indirizzo precedentemente menzionato) Iniziamo con un programma in C che prende in ingresso il file foto.ppm e andr`a a separare le componenti dell’immagine nei file matr R, matr G e matr B, oltre a careare un piccolo file contentente le dimensioni originali della foto. Ovviamente ciascuna componente non verr`a salvata come matrice, ma come sequenza di numeri: recupereremo la forma originale con il reshape() di Matlab. /* Programma che estrae dal file .ppm le componenti RGB separatamente. */ #include int main() { FILE *f1, *f2, *f3, *f4, *f5; int r, g, b; int alt, largh; f1 = fopen("foto.ppm", "r"); f2 = fopen("matr_R","w"); f3 = fopen("matr_G","w"); f4 = fopen("matr_B","w"); f5 = fopen("dimensione_foto","w"); fseek(f1,2*sizeof(char),SEEK_SET); 9

fscanf(f1,"%d %d", &largh, &alt); fprintf(f5,"%d %d", largh, alt); fclose(f5); fseek(f1,sizeof(int),SEEK_CUR); while (feof(f1) != 1) { fscanf(f1,"%d %d %d", &r, &g, &b); if (feof(f1) != 1) { fprintf(f2,"%d\n", r); fprintf(f3,"%d\n", g); fprintf(f4,"%d\n", b); } } fclose(f1); fclose(f2); fclose(f3); fclose(f4); return 0; } Abbiamo fatto uso di alcune operazioni su file: il comando fopen apre il file specificato in lettura se `e presente la stringa ’’r’’, in scirttura con la ’’w’’. I comandi fprintf e fscanf permettono di scrivere/leggere sul file indicato, mentre il comando fseek permette di spostare il cursore di lettura di un numero prefissato di bytes dalla posizione iniziale se specifico SEEK SET, dalla posizione attuale se specifico SEEK CUR, dalla posizione finale se specifico SEEK END. Salvo il file come crea matrici.c, non compilo in quanto lo far`o dopo direttamente da script. Scriviamo ora la function Matlab effettivamente adibita alla compressione, questa legger` a le componenti create col programma precedente, la dimensione dell’immagine e da tali dati effettuer`a le istruzioni elencate dal punto 1 al punto 3. % Programma che esegue compressione SVD su foto RGB. % k rappresenta la dimensione di ciascun sottoblocco. % m rappresenta il numero dei valori sing. da conservare. % z rappresenta il fattore di normalizzazione (cons. z = 20). function compress_svd(k,m,z) % Prima parte: adattamento componenti alla suddivisione in blocchi. disp(’Adatto componenti RGB...’) dim = load(’dimensione_foto’); largh = dim(1); alt = dim(2); foto_red = load(’matr_R’,’-ascii’); foto_green = load(’matr_G’,’-ascii’); foto_blue = load(’matr_B’,’-ascii’); foto_red = reshape(foto_red,largh,alt)’; foto_green = reshape(foto_green,largh,alt)’; foto_blue = reshape(foto_blue,largh,alt)’; 10

if (mod(largh,k) ~= 0) for i=1:k-mod(largh,k) foto_red(1:end,largh+i) = 128; foto_green(1:end,largh+i) = 128; foto_blue(1:end,largh+i) = 128; end end if (mod(alt,k) ~= 0) for i=1:k-mod(alt,k) foto_red(alt+i,1:end) = 128; foto_green(alt+i,1:end) = 128; foto_blue(alt+i,1:end) = 128; end end foto_red = foto_red(1:alt+k*(mod(alt,k) ~= 0)-mod(alt,k),1:largh+ + k*(mod(largh,k) ~= 0)-mod(largh,k)); foto_green = foto_green(1:alt+k*(mod(alt,k) ~= 0)-mod(alt,k),1:largh+ + k*(mod(largh,k) ~= 0)-mod(largh,k)); foto_blue = foto_blue(1:alt+k*(mod(alt,k) ~= 0)-mod(alt,k),1:largh+ + k*(mod(largh,k) ~= 0)-mod(largh,k)); sz = size(foto_red); dim(1) = sz(2); dim(2) = sz(1); % Seconda parte: estrazione svd e preparazione alla compressione. disp(’Comprimo immagine...’) matr_red = zeros((dim(1)/k)*(dim(2)/k),2*m); matr_green = zeros((dim(1)/k)*(dim(2)/k),2*m); matr_blue = zeros((dim(1)/k)*(dim(2)/k),2*m); cursore_matr = 1; for i=1:k:dim(1)-k+1 for j=1:k:dim(2)-k+1 % ROSSO [U,sigma,V] = svd(foto_red(j:j+k-1,i:i+k-1)); U = U(:,1:m); V = V(:,1:m); for p=1:m U(:,p) = U(:,p)*sqrt(sigma(p,p)); V(:,p) = V(:,p)*sqrt(sigma(p,p)); end matr_red(cursore_matr:cursore_matr+k-1,1:m) = floor(U*z); matr_red(cursore_matr:cursore_matr+k-1,m+1:2*m) = floor(V*z); % VERDE 11

[U,sigma,V] = svd(foto_green(j:j+k-1,i:i+k-1)); U = U(:,1:m); V = V(:,1:m); for p=1:m U(:,p) = U(:,p)*sqrt(sigma(p,p)); V(:,p) = V(:,p)*sqrt(sigma(p,p)); end matr_green(cursore_matr:cursore_matr+k-1,1:m) = floor(U*z); matr_green(cursore_matr:cursore_matr+k-1,m+1:2*m) = floor(V*z); % BLU [U,sigma,V] = svd(foto_blue(j:j+k-1,i:i+k-1)); U = U(:,1:m); V = V(:,1:m); for p=1:m U(:,p) = U(:,p)*sqrt(sigma(p,p)); V(:,p) = V(:,p)*sqrt(sigma(p,p)); end matr_blue(cursore_matr:cursore_matr+k-1,1:m) = floor(U*z); matr_blue(cursore_matr:cursore_matr+k-1,m+1:2*m) = floor(V*z); cursore_matr = cursore_matr+k; end end dim_matr = size(matr_red); save(’dim_matrice_compressa’,’dim_matr’,’-ascii’); file = fopen(’**compressa**’,’w’); fprintf(file,’%d\n’,matr_red); fprintf(file,’%d\n’,matr_green); fprintf(file,’%d\n’,matr_blue); fclose(file); end Vediamo alcuni punti nel programma da chiarire: le matrici vengono caricate da file con l’istruzione load, il cui parametro -ascii specifica a Matlab il fatto che le informazioni non sono codificate nel formato binario proprietario. La matrice di reshape va poi presa trasposta poich`e Matlab l’ha ricostruita per colonne, mentre l’immagine era invece codificata per righe. Infine, salviamo le dimensioni della matrice compressa M con il comando save, mentre salviamo la matrice stessa in forma vettoriale (in modo da poter risparmiare memoria evitando spazi vuoti) e per fare questo usiamo i comandi fopen fprintf fscanf in stile C. Abbiamo terminato la fase di compressione: non resta che scrivere lo script .sh in modo da far eseguire il tutto in automatico. #!/bin/bash echo "Dimensione dei sottoblocchi?" read NUM_QUAD 12

echo "Quanti valori singolari vuoi conservare?" read NUM_SV echo "Creo matrici per Octave/Matlab..." gcc -o crea_matrici crea_matrici.c ./crea_matrici echo "compress_svd(${NUM_QUAD},${NUM_SV},20)" | matlab echo "Elimino file non piu’ utili..." rm "matr_R" rm "matr_G" rm "matr_B" echo "Fatto!" Ora che la parte relativa alla compressione `e completa, vediamo come ricostruire un’immagine .ppm a partire dal file contenente i valori della matrice M . % Programma che decomprime l’immagine. % k rappresenta la dimensione dei sottoblocchi scelta in compressione. % z rappresenta il fattore di normalizzazione scelto in compressione. % t rappresente il valore di taglio dei vettori singolari. function decompress_svd(t,k,z) dim = load(’dim_matrice_compressa’,’-ascii’); %Carico elementi necessari. file = fopen(’**compressa**’,’r’); matr_red = fscanf(file,’%d’,dim(1)*dim(2)); matr_red = reshape(matr_red,dim(1),dim(2)); matr_green = fscanf(file,’%d’,dim(1)*dim(2)); matr_green = reshape(matr_green,dim(1),dim(2)); matr_blue = fscanf(file,’%d’,dim(1)*dim(2)); matr_blue = reshape(matr_blue,dim(1),dim(2)); fclose(file); dim_foto = load(’dimensione_foto’,’-ascii’); largh = dim_foto(1)+k*(mod(dim_foto(1),k) ~= 0)-mod(dim_foto(1),k); alt = dim_foto(2)+k*(mod(dim_foto(2),k) ~= 0)-mod(dim_foto(2),k); cursore_matr = 0; foto_red = zeros(alt,largh); %Inizializzo matrici dell’immagine. foto_green = zeros(alt,largh); foto_blue = zeros(alt,largh); blocchi_largh = largh/k; blocchi_alt = alt/k; for i=1:blocchi_largh for j=1:blocchi_alt for m=1:k for n=1:k somma = [0, 0, 0]; for o=1:t somma(1) = somma(1) + matr_red(m+cursore_matr,o)* matr_red(n+cursore_matr,t+o); 13

somma(2) = somma(2) + matr_green(m+cursore_matr,o)* matr_green(n+cursore_matr,t+o); somma(3) = somma(3) + matr_blue(m+cursore_matr,o)* matr_blue(n+cursore_matr,t+o); end val = floor(somma./(z^2)); val = max(0,val); val = min(255,val); foto_red(m+(j-1)*k,n+(i-1)*k) = val(1); foto_green(m+(j-1)*k,n+(i-1)*k) = val(2); foto_blue(m+(j-1)*k,n+(i-1)*k) = val(3); end end cursore_matr = cursore_matr+k; end end %Riporto la foto alle dimensioni originali. foto_red = foto_red(1:dim_foto(2),1:dim_foto(1)); %Calcolo trasposta, in modo da scrivere per riga. foto_red = foto_red’; foto_green = foto_green(1:dim_foto(2),1:dim_foto(1)); foto_green = foto_green’; foto_blue = foto_blue(1:dim_foto(2),1:dim_foto(1)); foto_blue = foto_blue’; %Salvo matrici. red = fopen(’matr_R’,’w’); green = fopen(’matr_G’,’w’); blue = fopen(’matr_B’,’w’); fprintf(red,’%d\n’,foto_red); fprintf(green,’%d\n’,foto_green); fprintf(blue,’%d\n’,foto_blue); fclose(red); fclose(green); fclose(blue); end Il codice sopra proposto non fa che eseguire i punti 4 e 5, ossia la ricostruzione della matrice approssimante blocco a blocco a partire dai valori singolari conservati e il recupero delle dimensioni originali contenute nel file dimensione foto creato in compressione. Alla fine di tale procedura resta da riformattare i dati in modo che l’immagine sia nuovamente riconoscibile come .ppm, in pratica `e sufficiente ricaricare i tre file in uscita dalla function e costruire un file unico introdotto dal valore di riconoscimento P3. Il tutto `e svolto da questo programma in C che salveremo come crea ppm.c. /* Programma che crea un’immagine .ppm a partire dalle matrici RGB. */ #include

14

int main() { FILE *f1,*f2,*f3,*f4,*f5; int r, g, b, alt, largh; f1 = fopen("matr_R","r"); f2 = fopen("matr_G", "r"); f3 = fopen("matr_B", "r"); f4 = fopen("dimensione_foto","r"); f5 = fopen("decompressa.ppm","w"); fscanf(f4,"%d %d", &largh, &alt); fclose(f4); fprintf(f5,"P3\n%d %d\n255\n", largh, alt); while (feof(f1) == 0) { fscanf(f1,"%d",&r); fscanf(f2,"%d",&g); fscanf(f3,"%d",&b); if (feof(f1) == 0) { fprintf(f5,"%d\n%d\n%d\n",r,g,b); } } fclose(f1); fclose(f2); fclose(f3); fclose(f5); return 0; } Siamo infine pronti a scrivere lo script di esecuzione e concludere cos`ı il codice relativo all’intero algoritmo. #!/bin/bash echo "Quanti valori singolari hai conservato in compressione?" read NUM_VAL echo "Di che dimensione sono i sottoblocchi?" read NUM_QUAD echo "Decomprimo l’immagine..." echo "decompress_svd(${NUM_VAL},${NUM_QUAD},20)" | matlab echo "Creo file .pgm..." gcc -o crea_ppm crea_ppm.c ./crea_ppm echo "Elimino files non piu’ utili..." rm "matr_R" rm "matr_G" rm "matr_B" echo "Fatto!" Passiamo ora alla sperimentazione vera e propria, considerando i risultati anche in relazione allo spazio occupato dal file di compressione. Per un buon rapporto qualit` a/memoria scegliamo p = 32: ci`o consente una procedura eseguita in tempi ragionevoli (ricordiamo che la singola SVD costa O(p3 )) e allo stesso tempo discreti risultati (ci `e sufficiente conservare solo 4 valori singolari per non notare 15

differenze alle dimensioni originali dell’immagine). Vogliamo essenzialmente contraddistinguere la variazione dell’algoritmo SVD da noi proposta con l’applicazione diretta sull’intera matrice immagine: dalle osservazioni emerge il fatto che suddividendo la foto in blocchi riesco a contenere il disturbo generale a scapito tuttavia della definizione dei bordi e di alcuni problemi sulle variazioni di tono. In questo modo la compressione ha un risultato minimo garantito, con una soglia minima di memoria occupata, mentre con il metodo SVD tradizionale si pu`o raggiungere il disturbo totale, non essendo pi` u in grado di distinguere i soggetti della figura. Raggiunto il grado ottimo di compressione per entrambi i metodi, non si nota un effettivo vantaggio in spazio di uno dei due: il lavoro pi` u preciso ottenuto con il nostro metodo lo paghiamo con una maggiore quantit` a di memoria. ` anzitutto necessario procurarsi un’imCome effettuiamo la sperimentazione? E magine loseless e questo `e un compito non cos`ı facile, visto che il 99% delle immagini sulla rete sono in realt` a gi`a compresse (proprio per gestire una migliore condivisione). Ci vengono in aiuto i produttori stessi di fotocamere che mettono a disposizione delle ”sample pictures” in RAW proprietario, la conversione in formato ASCII .ppm la effettuiamo con Gimp che con l’aiuto dell’add-on Ufraw riesce a gstire molti standard commerciali. Scegliamo in particolare la foto di esempio della Canon EOS 5D Mark III: l’immagine coniene notevoli variazioni di tono e buon dettaglio, nel formato .ppm occupa ben 49,7 MB! Rinominiamo il file come foto.ppm e collochiamolo nella stessa directory dove abbiamo salvato tutti i programmi precedenti. Eseguiamo varie compressioni: conservando 15 valori l’immagine compressa arriver` a a pesare 42,9 MB, con 10 valori 29,9 MB, con 4 valori 13,4 MB con 1 solo valore 4,6 MB. Siamo ovviamente abituati nel quotidiano a risultati molto migliori. Vediamo alcune immagini: poich`e a dimensioni originali non `e veramente possibile cogliere molte differenze (ovviamente sopra una certa soglia) pare pi` u opportuno osservare alcuni particolari opportunamente zoomati, scegliendo di focalizzarsi sulle sfumature chiare e frequenti ombreggiature della corolla di un fiore. Con questo concludiamo l’argomento.

16

Figura 4: l’immagine intera non compressa (.ppm)

Figura 5: particolari 15 sv, 4 sv, 3 sv, 1 sv.

17

La trasformata discreta del coseno (DCT) Sia v ∈ Rn , posso ovviamente esprimerlo come combinazione lineare dei vettori della base canonica {e1 , . . . , en }, dove   0  ..  .   0    ei =  1 (il valore 1 giace sull’i-esima riga) 0   .  ..  0 Ci` o significa che posso scrivere v come: v=

n X

αi ei .

i=1

Ovviamente αi = vi : ho enunciato questi principi estremamente banali con lo scopo di mostrare che la rappresentazione su base canonica, seppur fondamentale dal punto di vista operativo e rappresentativo, non fornisce informazioni utili al di fuori delle componenti stesse del vettore. Una base pi` u interessante a tal scopo `e quella costruita dai {c1 , . . . , cn }, definiti come:   √1 n

  c1 =  ...  √1 n



q

2 n

cos π(i−1) 2n



  q   3π(i−1) 2 cos   n 2n  , i = 2, . . . , n  ci =   ..   . q  π(2n−1)(i−1) 2 cos n 2n {c1 , . . . , cn } `e detta base del coseno; l’applicazione lineare Rn → Rn che associa un vettore v al vettore w dei coefficienti della combinazione su {c1 , . . . , cn } `e detta Trasformata discreta del coseno, in sigla DCT. Tra tutti i modelli formulati noi andremo sempre ad utilizzare la DCT-II che pu` o essere riassunta nella relazione che segue. wk = (DCTn (v))k = ζk

n X

vh cos

h=1

π(2h − 1)(k − 1) , k = 1, . . . , n 2n

dove,   √1 n ζk = q 2 

n

18

se k = 1 altrimenti

L’applicazione inversa, detta IDCT, associa il vettore dei coeffcienti w al vettore v ottenuto dalla combinazione lineare degli elementi di w sulla base {c1 , . . . , cn }. n X

vk = (IDCTn (w))k =

ζh wh cos

h=1

π(2k − 1)(h − 1) , k = 1, . . . , n 2n

Entrambe le procedure sono presenti di stock nella versione R2013a di Matlab: dct(v) calcola la trasformata diretta, mentre idct(w) l’inversa. Da segnalare il fatto che se effettuo la dct oppure l’idct su una matrice, Matlab effettuer`a l’applicazione sul vettore ottenuto concatenando le colonne della suddetta. Vediamo di estendere quanto detto finora sullo spazio vettoriale delle matrici Rn×n , poich`e `e qui che dovremmo applicare la trasformata per i nostri scopi di compressione. Se A ∈ Rn×n , la base canonica `e data da {E11 , . . . , E1n , E21 , . . . , E2n , . . . , En1 , . . . , Enn } dove,  0       Eij =  ...      0

... ..

0

     ..  , (l’1 giace in posizione (i, j)) .      0

. 0 1 0 ..



.

...

La base del coseno sulle matrici {C11 , . . . , C1n , E21 , . . . , C2n , . . . , Cn1 , . . . , Cnn } `e invece cos`ı strutturata:   π(2j − 1)(k − 1) π(2i − 1)(h − 1) cos Chk = ξh θk cos i = 1, . . . , n 2n 2n j = 1, . . . , n    √1  √1 se h = 1 se k = 1 n n dove, ξh = q 2 e θk = q 2 .   altrimenti altrimenti n n La DCT sar` a dunque un’applicazione lineare Rn×n → Rn×n e pu`o essere calcolata termine a termine con la seguente relazione: (B)ij = (DCTn×n (A))ij = ξi θj

n X n X

(A)hk cos

h=1 k=1

π(2k − 1)(j − 1) π(2h − 1)(i − 1) cos 2n 2n

Analogo ragionamento per la trasformata inversa. (A)ij = (IDCTn×n (B))ij =

n X n X

ξh θk (B)hk cos

h=1 k=1

19

π(2i − 1)(h − 1) π(2j − 1)(k − 1) cos 2n 2n

Come possiamo utilizzare la DCT matriciale e la sua inversa? Sfruttando la function dct possiamo fabbricarci la nostra my dct2 conoscendo la trasformata per colonne della matrice e della rispettiva trasposta. Ci risparmiamo lo sforzo, visto che Matlab integra le funzioni dct2 e idct2 che fanno tutto il lavoro ”sporco” per noi! Matlab inoltre, mediante algorimti efficienti, riesce ad eseguire il tutto in O(n log n).

Figura 6: la base del coseno in R8×8 In quest’ultima parte del paragrafo parliamo del motivo fondamentale per cui si `e scelto di lavorare sulla DCT e del legame con l’algoritmo JPEG. La trasformata del coseno permette di valutare l’immagine secondo determinate variazioni di tono (come mostrato in Figura 6), conoscendo per ogni tipo di variazione il coefficiente di intensit`a. Le variazioni di tono pi` u repentine sono dette ”frequenze alte”, metre quelle pi` u graduali sono le ”frequenze basse”: l’occhio umano tende a focalizzarsi sulle alte frequenze trascurando le basse. Considerando anche il fatto che in un’immagine le variazioni graduali avranno inoltre coefficiente pi` u alto, `e evidente che la DCT sar`a il nucleo imprescindibile di un algoritmo di compressione JPEG. Facciamo un esperimento: prendiamo una foto .ppm, acquisiamola con Matlab e calcoliamone la DCT. Notiamo che le basse frequenze si addenseranno in alto a sinistra (vedi Figura 8), baster`a eliminare le altre e avremo la nostra compressione; di questo ne disuteremo con abbondanza nel successivo paragrafo.

20

Figura 7: Immagine scelta come campione.

Figura 8: un’imagesc della DCT eseguita sulla matrice dei rossi. I valori pi` u elevati in alto a sinistra.

21

Secondo approccio: l’algoritmo JPEG. Abbiamo gi` a introdotto le conoscenze necessarie, entriamo subito nel vivo enunciando le fasi della compressione JPEG: 1. Separazione dell’immagine in blocchi p×p, di default p = 8, e calcolo della DCT su ogni blocco. 2. Quantizzazione, ossia la tecnica di cernita nei confronti dei coefficenti della DCT corrispondenti alle alte frequenze. 3. Codifica, ovvero la costruzione effettiva del file .jpg che noi non tratteremo poich`e lontana dai nostri fini di sperimentazione. Prima di procedere con la spiegazione dettagliata dei punti precedenti, alcune riflessioni sono d’obbligo. • Abbiamo abbondantemente preannunciato che il JPEG si esegue sullo spazio colore di lumiocrominanza e non su RGB: questa oltre ad essere una standardizzazione permette di avere risultati migliori dal punto di vista del colore e delle sfumature, in quanto la compressione `e eseguita autonomamente sulla componente luminosa rispetto a quella cromatica. L’YCbCr consente un miglior fattore di risparmio memoria a parit`a di qualit` a visiva, a scapito tuttavia di una computazione meno rapida (sar`a necessario infatti scrivere del codice di conversione prima e dopo la fase di compressione). Una curiosit`a derivante dall’uso di tale spazio colore consiste nel differente degrado dell’immagine derivante da eccessiva severit` a di compressione: anzich`e distrurbare la sagoma degli oggetti, come nel RGB, l’algoritmo andr`a a colpire le variazioni, rendendo visibili delle vere e proprie ”curve di livello”. • Cerchiamo di motivare la scelta di p = 8. Una suddivisione pi` u fine porterebbe ad un risultato pi` u accurato, ma ad un tempo di esecuzione non propriamente sostenibile; allo stesso tempo una scelt`a pi` u larga non sarebbe esente da difetti visivi sull’immagine compressa. Scegliendo p = 32 notiamo che il computer riesce a portare a termine il tutto molto pi` u velocemente, ma con facilit`a insorgono problemi. Supponiamo che un blocco includa un tratto con alte variazioni (ad es. parte di un petalo del fiore) e il resto a variazione pi` u lieve (ad es. lo sfondo), mentre il blocco adiacente sia esclusivamente a bassa variazione: sar` a visibile lo ”sfarfallio” proveniente dalle alte frequenze in contrasto con una tonalit` a uniforme del secondo blocco. p = 8 rappresenta ancora una volta un buon compromesso. • Dal punto di vista matematico sappiamo che l’approccio SVD `e una migliore approssimazione, ma come vedremo il JPEG consente un risparmio di spazio nettamente maggiore con eguale qualit`a. Quest’algoritmo non si basa infatti sul ricostruire l’immagine di partenza il pi` u fedelmente possibile, ma sull’eliminazione di quelle componenti dell’immagine inutili perch`e non percettibili. L’esperienza vince dunque sulla teoria.

22

Descrivendo pi` u precisamente il nostro elenco numerato, quali coefficenti della DCT andremo ad eliminare e quali conserveremo? Un’idea ingenua potrebbe essere l’eliminazione di tutti quei valori il cui modulo `e al di sotto di una data soglia: avremmo risultati gi`a notevoli in quanto, come abbiamo visto in Figura 8, le basse frequenze corrispondono ai valori pi` u elevati. Una compressione svolta in questo modo non sarebbe il massimo dell’accuratezza, ma il vero problema `e che non sapremmo quantifcare l’effettivo grado di riduzione, dovremmo stabilire una soglia di taglio specifica in base alla foto. La soluzione ortodossa `e data dal processo di quantizzazione: il JPEG, una volta ottenuta la DCT del blocco, esegue una divisione delle componenti termine a termine con delle matrici fornite opportunamente riscalate. I valori decimali vengono arrotondati e, come vedremo, molti di essi saranno ridotti a 0. Le matrici di quantizzazione dello standard JPEG sono le seguenti, la prima per la componente di lumia, mentre la seconda per le due componenti di crominanza:   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    QY =   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

QCb,Cr

 17 18  24  47 = 99  99  99 99

18 21 26 66 99 99 99 99

24 26 56 99 99 99 99 99

47 66 99 99 99 99 99 99

99 99 99 99 99 99 99 99

99 99 99 99 99 99 99 99

99 99 99 99 99 99 99 99

 99 99  99  99  99  99  99 99

Non esiste un metodo preciso per determinare le suddette, se non per procedure euristiche: la scelta delle matrici di quantizzazione non `e unica e tutt’oggi si lavora per crearne sempre di migliori. Ciascun dispositivo o software in grado di eseguire il JPEG ne ha di proprie, possiamo trovarne alcune all’indirizzo http://www.impulseadventure.com/photo/jpeg-quantization.html. Si osserva che, in conformit` a alla Figura 8, le matrici di quantizzazione hanno tutte la medesima struttura: i valori in alto a sinistra sono pi` u piccoli, in modo che la cancellazione agisca maggiormente sulle altre aree del blocco. Come anticipato prima, non lavoreremo sulle matrici QY e QCb,Cr direttamente: fornito un intero l ∈ [1, 100] (1 per la massima compressione possibile), riscaleremo le matrici di quantizzazione secondo le seguenti regole, ottenendo Q0Y e Q0Cb,Cr . ( 5000 se l < 50 l µ= 200 − 2l altrimenti Q0ij = b

µQij + 50 c 100

23

A questo punto l’algoritmo JPEG procede con la memorizzazione dei coefficenti secondo opportune codifiche esplicitando i dati in esadecimale: tutto ci`o `e svolto secondo un’opportuno standard (detto ”codifica entropica”) che consente di far aprire il file compresso direttamente ai visualizzatori di immagine, senza implementare decompressioni. Noi trascureremo questa fase, ricorrendo alla gestione abbastanza vantaggiosa delle matrici sparse su Matlab. Andiamo a descrivere il nostro algoritmo: analogamente a quanto svolto per la SVD, da ciascuna componente A dell’immagine costruiamo la A˜ in modo che sia p × p-divisibile. Eseguiamo l’esplorazione blocco a blocco, procedendo per colonne: per ciascun blocco B(i,j) calcoliamo K(i,j) = q[DCTn×n (B(i,j) )], q[.] rappresenta l’operazione di quantizzazione prima descritta. Memorizziamo la mn matrice M ∈ R p ×p cos`ı strutturata:   K(1,1)  K(2,1)      ...      K M =  (1,2)    K(2,2)      ..   . K( mp , np ) Vista la larghissima presenza di zeri, useremo l’istruzione sparse(M) di Matlab in modo da salvare solo gli elementi non nulli con la loro effettiva posizione. Il fatto che M abbia notevoli dimensioni in lunghezza consente all’istruzione sparse di far risparmiare un’enorme quantit`a di memoria, alla stregua della codifica standard. Poich`e non ci `e possbile leggere il file di compressione come immagine, `e necessario implementare la decompressione: per ciascun elemento K(i,j) andremo a ricavarci B 0 (i,j) = q −1 [IDCTn×n (K(i,j) )], dove q −1 [.] `e l’operazione di dequantizzazione, ossia la moltiplicazione del blocco termine a termine per gli elementi della matrice di quantizzazione. Ricostruisco A˜0 = (B 0 (i,j) )ij e riporto l’immagine alle dimensioni originali avendo infine A0 . Osserviamo il codice in ogni sua parte: per separare le componenti RGB usiamo lo stesso file crea matrici.c scritto per la SVD, mentre implementiamo una function per la conversione nello spazio colore YCbCr. function RGB_to_YCbCr() % Carico matrici delle componenti RGB. disp(’Carico componenti RGB...’); matr_R = load(’matr_R’,’-ascii’)-128; matr_G = load(’matr_G’,’-ascii’)-128; matr_B = load(’matr_B’,’-ascii’)-128; % Carico dimensioni della foto. dim = load(’dimensione_foto’,’-ascii’); % Inizializzo matrici YCbCr. matr_Y = zeros(dim(2),dim(1)); 24

matr_Cb = zeros(dim(2),dim(1)); matr_Cr = zeros(dim(2),dim(1)); % Matrice di trasformazione. A = [0.299 0.587 0.114; -0.1687 -0.3313 0.5; 0.5 -0.4187 -0.0813]; % Eseguo appl. lineare. disp(’Eseguo conversione a YCbCr...’); for i=1:dim(1)*dim(2) vett = [matr_R(i); matr_G(i); matr_B(i)]; vett = A*vett; matr_Y(i) = vett(1); matr_Cb(i) = vett(2); matr_Cr(i) = vett(3); end % Salvo matrici YCbCr. disp(’Salvo dati...’); Y = fopen(’matr_Y’,’w’); Cb = fopen(’matr_Cb’,’w’); Cr = fopen(’matr_Cr’,’w’); fprintf(Y,’%d\n’,matr_Y); fprintf(Cb,’%d\n’,matr_Cb); fprintf(Cr,’%d\n’,matr_Cr); fclose(Y); fclose(Cb); fclose(Cr); end In questa procedura, come gi` a spiegato all’inizio, abbiamo shiftato tutte le componenti RGB di -128, dopodich`e per ogni tripla si `e eseguita l’applicazione lineare di passaggio alla lumiocrominanza. Abbiamo salvato le componenti in forma vettoriale nei file matr Y, matr Cb e matr Cr. Dobbiamo ora creare un file contenente le nostre matrici di quantizzazione che chiameremo matrici quantizzazione.mat, per fare ci`o usiamo il formato binario proprietario di Matlab. >>> >>> >>> >>>

quant_lum = %Riporto la matrice Q_Y; quant_crom = %Riporto la matrice Q_Cb,Cr; save(’matrici_quantizzazione.mat’,’quant_lum’); save(’matrici_quantizzazione.mat’,’quant_crom’,’-append’);

Abbiamo tutto il necessario per procedere alla compressione. % Programma che esegue una compressione DCT sulle componenti YCbCr, % mediante quantizzazione. % k rappresenta la dimensione di ciacun sottoblocco in cui viene divisa la foto. % scal ` e il fattore di scalatura della matrice di quantizzazione, compreso tra 1 e 100. 25

function compress_dct(k,scal) % Prima parte: scalatura delle matrici di quantizzazione. disp(’Scalatura matrici di quantizzazione...’); load(’matrici_quantizzazione.mat’); if (scal < 50) quant_lum = ceil((quant_lum.*(5000/scal)+50*ones(8))./100); quant_crom = ceil((quant_crom.*(5000/scal)+50*ones(8))./100); else quant_lum = ceil((quant_lum.*(200-scal*2)+50*ones(8))./100); quant_crom = ceil((quant_crom.*(200-scal*2)+50*ones(8))./100); end save(’quant_scalata.mat’,’quant_lum’); % Salvo scalatura. save(’quant_scalata.mat’,’quant_crom’,’-append’); % Seconda parte: adattamento immagine alla suddivisione in blocchi. disp(’Adatto componenti...’) dim = load(’dimensione_foto’,’-ascii’); largh = dim(1); alt = dim(2); foto_Y = load(’matr_Y’,’-ascii’); foto_Cb = load(’matr_Cb’,’-ascii’); foto_Cr = load(’matr_Cr’,’-ascii’); foto_Y = reshape(foto_Y,largh,alt)’; foto_Cb = reshape(foto_Cb,largh,alt)’; foto_Cr = reshape(foto_Cr,largh,alt)’; if (mod(largh,k) ~= 0) % Adattamento larghezza foto. for i=1:k-mod(largh,k) foto_Y(1:end,largh+i) = 128; foto_Cb(1:end,largh+i) = 128; foto_Cr(1:end,largh+i) = 128; end end if (mod(alt,k) ~= 0) % Adattamento altezza foto. for i=1:k-mod(alt,k) foto_Y(alt+i,1:end) = 128; foto_Cb(alt+i,1:end) = 128; foto_Cr(alt+i,1:end) = 128; end end foto_Y = foto_Y(1:alt+k*(mod(alt,k) ~= 0)-mod(alt,k),1:largh+ +k*(mod(largh,k) ~= 0)-mod(largh,k)); foto_Cb = foto_Cb(1:alt+k*(mod(alt,k) ~= 0)-mod(alt,k),1:largh+ +k*(mod(largh,k) ~= 0)-mod(largh,k)); 26

foto_Cr = foto_Cr(1:alt+k*(mod(alt,k) ~= 0)-mod(alt,k),1:largh+ +k*(mod(largh,k) ~= 0)-mod(largh,k)); sz = size(foto_Y); dim(1) = sz(2); dim(2) = sz(1); % Terza parte: calcolo della DCT e compressione. disp(’Calcolo DCT e comprimo...’); matr_Y = zeros((dim(1)/k)*(dim(2)/k),k); matr_Cb = zeros((dim(1)/k)*(dim(2)/k),k); matr_Cr = zeros((dim(1)/k)*(dim(2)/k),k); cursore_matr = 1; for i=1:k:dim(1)-k+1 for j=1:k:dim(2)-k+1 trasf_Y = dct2(foto_Y(j:j+k-1,i:i+k-1)); % Calcolo trasformate. trasf_Cb = dct2(foto_Cb(j:j+k-1,i:i+k-1)); trasf_Cr = dct2(foto_Cr(j:j+k-1,i:i+k-1)); trasf_Y = round(trasf_Y./quant_lum); % Eseguo quantizzazione. trasf_Cb = round(trasf_Cb./quant_crom); trasf_Cr = round(trasf_Cr./quant_crom); matr_Y(cursore_matr:cursore_matr+k-1,:) = trasf_Y; matr_Cb(cursore_matr:cursore_matr+k-1,:) = trasf_Cb; matr_Cr(cursore_matr:cursore_matr+k-1,:) = trasf_Cr; cursore_matr = cursore_matr+k; end end disp(’Scrivo i dati...’); matr_Y = sparse(matr_Y); % Acquisizione come matrice sparsa. matr_Cb = sparse(matr_Cb); matr_Cr = sparse(matr_Cr); save(’compressa.mat’,’matr_Y’); save(’compressa.mat’,’matr_Cb’,’-append’); save(’compressa.mat’,’matr_Cr’,’-append’); end Nella prima parte eseguiamo la scalatura delle matrici memorizzate come spiegato in precedenza; unico accorgimento da segnalare `e l’uso di ceil, ossia l’arrotondamento all’intero superiore, in modo da non avere zeri in quantizzazione qualora si scelga un fattore di compressione pari a 100. La seconda parte, il ridimensionamento, `e perfettamente analoga a quanto svolto per l’SVD. Nella terza parte abbiamo eseguito per ogni unit`a di tutte le componenti la DCT, la quantizzazione e abbiamo salvato infine il tutto nel file compressa.mat. Da notare l’uso di -append per inserire pi` u matrici nello stesso file binario. La fase di compressione `e completata: scriviamo lo scipt .sh. #!/bin/bash echo "Livello di compressione? (1-100)" 27

read SCALA echo "Creo matrice per Octave/Matlab..." gcc -o crea_matrici crea_matrici.c ./crea_matrici echo "RGB_to_YCbCr();compress_dct(8,${SCALA})" | matlab echo "Elimino file non piu’ utili..." rm "matr_R" rm "matr_G" rm "matr_B" rm "matr_Y" rm "matr_Cb" rm "matr_Cr" echo "Fatto!" Lo script di compressione produce i file compressa.mat e dimensione foto a partire da foto.pgm: procediamo ora a creare il codice di decompressione. % Programma che esegue la decompressione calcolando per % ogni blocco la trasformata del coseno inversa a % partire dalla matrice sparsa dei coefficienti conservati. % k rappresenta la suddivisione dei sottoblocchi (k=8 per il JPEG). function decompress_dct(k) %Carico elementi necessari. load(’compressa.mat’); load(’quant_scalata.mat’); dim_foto = load(’dimensione_foto’,’-ascii’); %Recupero della forma completa. matr_Y = full(matr_Y); matr_Cb = full(matr_Cb); matr_Cr = full(matr_Cr);

largh = dim_foto(1)+k*(mod(dim_foto(1),k) ~= 0)-mod(dim_foto(1),k); alt = dim_foto(2)+k*(mod(dim_foto(2),k) ~= 0)-mod(dim_foto(2),k); cursore_matr = 1; %Inizializzo matrice dell’immagine. foto_Y = zeros(alt,largh); foto_Cb = zeros(alt,largh); foto_Cr = zeros(alt,largh); blocchi_largh = largh/k; blocchi_alt = alt/k; disp(’Decomprimo la foto...’); for i=1:blocchi_largh for j=1:blocchi_alt %Calcolo l’antitrasformata del coseno. lum = matr_Y(cursore_matr:cursore_matr+k-1,:); 28

lum = lum.*quant_lum; crom1 = matr_Cb(cursore_matr:cursore_matr+k-1,:); crom1 = crom1.*quant_crom; crom2 = matr_Cr(cursore_matr:cursore_matr+k-1,:); crom2 = crom2.*quant_crom; antitrasf_Y = round(idct2(lum)); antitrasf_Cb = round(idct2(crom1)); antitrasf_Cr = round(idct2(crom2)); foto_Y((j-1)*k+1:j*k,(i-1)*k+1:i*k) = antitrasf_Y; foto_Cb((j-1)*k+1:j*k,(i-1)*k+1:i*k) = antitrasf_Cb; foto_Cr((j-1)*k+1:j*k,(i-1)*k+1:i*k) = antitrasf_Cr; cursore_matr = cursore_matr+k; end end %Riporto la foto alle dimensioni originali. foto_Y = foto_Y(1:dim_foto(2),1:dim_foto(1)); foto_Cb = foto_Cb(1:dim_foto(2),1:dim_foto(1)); foto_Cr = foto_Cr(1:dim_foto(2),1:dim_foto(1)); %Calcolo trasposta, in modo da scrivere per riga. foto_Y = foto_Y’; foto_Cb = foto_Cb’; foto_Cr = foto_Cr’; %Salvo matrice. Y = fopen(’matr_Y’,’w’); Cb = fopen(’matr_Cb’,’w’); Cr = fopen(’matr_Cr’,’w’); fprintf(Y,’%d\n’,foto_Y); fprintf(Cb,’%d\n’,foto_Cb); fprintf(Cr,’%d\n’,foto_Cr); fclose(Y); fclose(Cb); fclose(Cr); end Questa function esegue, iterando su ogni blocco, la dequantizzazione e il calcolo dell’IDCT. In modo da avere infine solo valori interi uso l’arrotondamento svolto da round. Salvo i dati in componenti separate sui file matr Y, matr Cb e matr Cr: poich`e il file .ppm non pu` o leggere da lumiocrominanza devo rieseguire una conversione in RGB, faccio ci` o con il seguente. function YCbCr_to_RGB() % Carico matrici delle componenti YCbCr. disp(’Carico componenti YCbCr...’); matr_Y = load(’matr_Y’,’-ascii’); matr_Cb = load(’matr_Cb’,’-ascii’); matr_Cr = load(’matr_Cr’,’-ascii’); % Carico dimensioni della foto. 29

dim = load(’dimensione_foto’,’-ascii’); % Inizializzo matrici RGB. matr_R = zeros(dim(2),dim(1)); matr_G = zeros(dim(2),dim(1)); matr_B = zeros(dim(2),dim(1)); % Matrice di trasformazione. A = [1 0 1.402; 1 -0.344 -0.714; 1 1.772 0]; % Eseguo appl. lineare. disp(’Eseguo conversione a RGB...’); for i=1:dim(1)*dim(2) vett = [matr_Y(i); matr_Cb(i); matr_Cr(i)]; vett = A*vett; matr_R(i) = round(vett(1)+128); matr_G(i) = round(vett(2)+128); matr_B(i) = round(vett(3)+128); end % Controllo che i dati siano nel range [0,255]. matr_R = max(0,matr_R); matr_R = min(255,matr_R); matr_G = max(0,matr_G); matr_G = min(255,matr_G); matr_B = max(0,matr_B); matr_B = min(255,matr_B); % Salvo matrici RGB. disp(’Salvo dati...’); R = fopen(’matr_R’,’w’); G = fopen(’matr_G’,’w’); B = fopen(’matr_B’,’w’); fprintf(R,’%d\n’,matr_R); fprintf(G,’%d\n’,matr_G); fprintf(B,’%d\n’,matr_B); fclose(R); fclose(G); fclose(B); end Resta solo da ricostruire il .ppm, basta riutilizzare il programma C che abbiamo gi` a salvato come crea ppm.c; per quanto riguarda lo script di esecuzione #!/bin/bash echo "decompress_dct(8);YCbCr_to_RGB()" | matlab echo "Creo file .pgm..." 30

gcc -o crea_ppm crea_ppm.c ./crea_ppm echo "Elimino files non piu’ utili..." rm "matr_Y" rm "matr_Cb" rm "matr_Cr" rm "matr_R" rm "matr_G" rm "matr_B" echo "Fatto!" Siamo pronti a testare il nostro JPEG: in modo da avere un confronto pulito con l’algoritmo SVD precedente useremo la medeisima foto (49,7 MB). Notiamo che dal punto di vista del risparmio di spazio non ci sono termini di paragone: conservando 4 valori singolari avevamo il rapporto ottimale compressione/qualit` a con 13,4 MB, col JPEG ci`o accade impostando un fattore di scalatura 80 e solo 888 KB! Il risultato `e veramente eccellente, visto che si tratta di un’immagine ad alta risoluzione, 2808x1784 pixel. Vediamo invece come il JPEG deteriora l’immagine con l’aumentare dei coefficenti portati a 0: si osserva una perdita generale di qualit`a, il fatto che in maniera sparsa si iniziano a notare quadratini e macchie di colore, nei casi pi` u estremi si ha pure una perdita di naturalezza cromatica. Nell’algoritmo SVD notavamo invece distorsioni nei pressi dei bordi, ma mantenimento sulle aree uniformi: il JPEG ha un comportamento pi` u ”indifferente” rispetto alla natura dell’immagine in oggetto. ` sorprendente il fatto che con la conservazione di pochi coefficienti del coseE no, e quindi matrici di compressione fortemente sparse, si riesca a riscostruire un immagine altamente fedele all’originale: la DCT infatti ”cattura” e separa i dettagli importanti per la visione da quelli trascurabili e pertanto eliminabili. Inoltre, stabilita la dimensione 8x8 di suddivisione, che per pi` u ragioni `e ritenuta adatta, i coefficienti rappresentanti le componenti essenziali sono veramente pochi e ci` o non sarebbe effettivamente fattibile con suddivisioni di grandezza magggiore; si andrebbe a perdere in nitidezza. Nella nostra sperimentazione abbiamo ottenuto una dimensione di 535 KB con fattore di sclatura 60, 374 KB con fattore 30, 253 KB con fattore 12, 191 KB con fattore 7, 115 KB con fattore 3: da segnalare la scarsa qualit`a degli ultimi tre tentativi.

31

Figura 9: l’immagine intera non compressa (.ppm)

Figura 10: particolari con fattori di scalatura 60, 30, 12, 7.

32

Un esperimento: variazione della scalatura di quantizzazione. L’implementazione svolta nel paragrafo precedente consente di effettuare la compressione di un immagine a partire da un valore numerico nel range [1, 100]: per valori vicini a 1 otteniamo massimo risparmio di spazio con peggiore qualit`a, man mano che il valore aumenta crescono resa e spazio. Proviamo ora a far variare la matrice di quantizzazione Q in corso d’opera, in modo da testare pi` u m×n gradi di compressione sulla stessa foto: pi` u precsiamente, se A ∈ R sta  bilisco un vettore l1 . . . lr di livelli di compressione. Andr`o a scalare Q mn mn rispetto al fattore l1 nei primi 2 blocchi, rispetto a l2 nei successivi 2 e p r p r cos`ı via. L’algortimo che presenteremo ha tuttavia delle problematiche qualora r - mn, poich`e avremmo problemi di divisione non intera: supponiamo che la dimensione del nostro vettore in input divida il prodotto di altezza e larghezza foto. Ci `e necessario modificare solo le function di compressione e decompressione, il resto sar` a prettamente identico a quanto visto nel paragrafo precedente. % Programma che esegue la compressione dell’immagine mediante % DCT facendo variare la scalatura della matrice di quantizzazione. % v ` e il vettore contenente i fattori di compressione al singolo tratto. % Assumiamo length(v) divisibile per il numero di blocchi, per non aver % problemi di fraz. non intero. % k ` e la dimensione dei sottoblocchi (k=8 per il JPEG). function compress_tratti(v,k) % Carico file necessari. load(’matrici_quantizzazione.mat’); dim = load(’dimensione_foto’,’-ascii’); largh = dim(1); alt = dim(2); foto_Y = load(’matr_Y’,’-ascii’); foto_Cb = load(’matr_Cb’,’-ascii’); foto_Cr = load(’matr_Cr’,’-ascii’); foto_Y = reshape(foto_Y,largh,alt)’; foto_Cb = reshape(foto_Cb,largh,alt)’; foto_Cr = reshape(foto_Cr,largh,alt)’; % Adattamento immagine alla suddivisione in blocchi. if (mod(largh,k) ~= 0) % Adattamento larghezza foto. for i=1:k-mod(largh,k) foto_Y(1:end,largh+i) = 128; foto_Cb(1:end,largh+i) = 128; foto_Cr(1:end,largh+i) = 128; end end if (mod(alt,k) ~= 0) % Adattamento altezza foto. 33

for i=1:k-mod(alt,k) foto_Y(alt+i,1:end) = 128; foto_Cb(alt+i,1:end) = 128; foto_Cr(alt+i,1:end) = 128; end end foto_Y = foto_Y(1:alt+k*(mod(alt,k) ~= 0)-mod(alt,k),1:largh+ +k*(mod(largh,k) ~= 0)-mod(largh,k)); foto_Cb = foto_Cb(1:alt+k*(mod(alt,k) ~= 0)-mod(alt,k),1:largh+ +k*(mod(largh,k) ~= 0)-mod(largh,k)); foto_Cr = foto_Cr(1:alt+k*(mod(alt,k) ~= 0)-mod(alt,k),1:largh+ +k*(mod(largh,k) ~= 0)-mod(largh,k)); sz = size(foto_Y); dim(1) = sz(2); dim(2) = sz(1);

% Calcolo il numero di blocchi per tratto di compressione. num = ((dim(1)/k)*(dim(2)/k))/length(v); % Calcolo della DCT e compressione. disp(’Calcolo DCT e comprimo...’); matr_Y = zeros((dim(1)/k)*(dim(2)/k),k); matr_Cb = zeros((dim(1)/k)*(dim(2)/k),k); matr_Cr = zeros((dim(1)/k)*(dim(2)/k),k); cursore_matr = 1; prossima_riscalatura = 1; riscala = 1; % Var. logica che indica necessit` a % di cambiare scala. curr = 1; % Indicatore del tratto corrente. num_blocchi = 0; % Contatore di blocchi compressi. for i=1:k:dim(1)-k+1 for j=1:k:dim(2)-k+1 if (riscala == 1) % Riscalatura della matrice di quant. if (v(curr) < 50) quant_lum_new = ceil((quant_lum.*(5000/v(curr))+50*ones(8))./100); quant_crom_new = ceil((quant_crom.*(5000/v(curr))+50*ones(8))./100); else quant_lum_new = ceil((quant_lum.*(200-v(curr)*2)+50*ones(8))./100); quant_crom_new = ceil((quant_crom.*(200-v(curr)*2)+50*ones(8))./100); end riscala = 0; curr = curr+1; prossima_riscalatura = prossima_riscalatura+num; end trasf_Y = dct2(foto_Y(j:j+k-1,i:i+k-1)); trasf_Cb = dct2(foto_Cb(j:j+k-1,i:i+k-1)); 34

trasf_Cr = dct2(foto_Cr(j:j+k-1,i:i+k-1)); trasf_Y = round(trasf_Y./quant_lum_new); trasf_Cb = round(trasf_Cb./quant_crom_new); trasf_Cr = round(trasf_Cr./quant_crom_new); matr_Y(cursore_matr:cursore_matr+k-1,:) = trasf_Y; matr_Cb(cursore_matr:cursore_matr+k-1,:) = trasf_Cb; matr_Cr(cursore_matr:cursore_matr+k-1,:) = trasf_Cr; cursore_matr = cursore_matr+k; num_blocchi = num_blocchi+1; if (num_blocchi == prossima_riscalatura-1) riscala = 1; end end end disp(’Scrivo i dati...’); matr_Y = sparse(matr_Y); matr_Cb = sparse(matr_Cb); matr_Cr = sparse(matr_Cr); save(’compressa.mat’,’matr_Y’); save(’compressa.mat’,’matr_Cb’,’-append’); save(’compressa.mat’,’matr_Cr’,’-append’); end Metodicamente non vi sono differenze con l’algoritmo JPEG originale; `e utilizzata la variabile logica riscala che segnala la necessita di modificare la matrice di quantizzazione una volta raggiunto un numero di blocchi esaminati pari a prossima riscalatura, quest’ultima variabile `e prontamente aggiornata ogni volta. Identici accorgimenti sono utilizzati nella function di decompressione: non sono dunque necessarie ulteriori spiegazioni. % Programma che esegue la decompressione della foto. % v ` e il vettore contenente i fattori di compressione al singolo tratto. % Assumiamo length(v) divisibile per il numero di blocchi, per non aver % problemi di fraz. non intero. % k ` e la dimensione dei sottoblocchi (k=8 per il JPEG). function decompress_tratti(v,k) load(’compressa.mat’); %Carico elementi necessari. load(’matrici_quantizzazione.mat’); dim_foto = load(’dimensione_foto’,’-ascii’); matr_Y = full(matr_Y); %Recupero della forma completa. matr_Cb = full(matr_Cb); matr_Cr = full(matr_Cr);

largh = dim_foto(1)+k*(mod(dim_foto(1),k) ~= 0)-mod(dim_foto(1),k); alt = dim_foto(2)+k*(mod(dim_foto(2),k) ~= 0)-mod(dim_foto(2),k); cursore_matr = 1; 35

prossima_riscalatura = 1; riscala = 1; curr = 1; num_blocchi = 0; foto_Y = zeros(alt,largh); %Inizializzo matrice dell’immagine. foto_Cb = zeros(alt,largh); foto_Cr = zeros(alt,largh); blocchi_largh = largh/k; blocchi_alt = alt/k; num = (blocchi_largh*blocchi_alt)/length(v); disp(’Decomprimo la foto...’); for i=1:blocchi_largh for j=1:blocchi_alt %Calcolo l’antitrasformata del coseno. if (riscala == 1) % Riscalatura della matrice di quant. if (v(curr) < 50) quant_lum_new = ceil((quant_lum.*(5000/v(curr))+50*ones(8))./100); quant_crom_new = ceil((quant_crom.*(5000/v(curr))+50*ones(8))./100); else quant_lum_new = ceil((quant_lum.*(200-v(curr)*2)+50*ones(8))./100); quant_crom_new = ceil((quant_crom.*(200-v(curr)*2)+50*ones(8))./100); end riscala = 0; curr = curr+1; prossima_riscalatura = prossima_riscalatura+num; end lum = matr_Y(cursore_matr:cursore_matr+k-1,:); lum = lum.*quant_lum_new; crom1 = matr_Cb(cursore_matr:cursore_matr+k-1,:); crom1 = crom1.*quant_crom_new; crom2 = matr_Cr(cursore_matr:cursore_matr+k-1,:); crom2 = crom2.*quant_crom_new; antitrasf_Y = round(idct2(lum)); antitrasf_Cb = round(idct2(crom1)); antitrasf_Cr = round(idct2(crom2)); foto_Y((j-1)*k+1:j*k,(i-1)*k+1:i*k) = antitrasf_Y; foto_Cb((j-1)*k+1:j*k,(i-1)*k+1:i*k) = antitrasf_Cb; foto_Cr((j-1)*k+1:j*k,(i-1)*k+1:i*k) = antitrasf_Cr; cursore_matr = cursore_matr+k; num_blocchi = num_blocchi+1; if (num_blocchi == prossima_riscalatura-1) riscala = 1; end end end foto_Y = foto_Y(1:dim_foto(2),1:dim_foto(1)); 36

foto_Cb = foto_Cb(1:dim_foto(2),1:dim_foto(1)); foto_Cr = foto_Cr(1:dim_foto(2),1:dim_foto(1)); foto_Y = foto_Y’; foto_Cb = foto_Cb’; foto_Cr = foto_Cr’; Y = fopen(’matr_Y’,’w’); Cb = fopen(’matr_Cb’,’w’); Cr = fopen(’matr_Cr’,’w’); fprintf(Y,’%d\n’,foto_Y); fprintf(Cb,’%d\n’,foto_Cb); fprintf(Cr,’%d\n’,foto_Cr); fclose(Y); fclose(Cb); fclose(Cr); end Proviamo quanto abbiamo implementato su un’immagine piuttosto famosa: trattasi di una modella svedese, Lena Soderberg, la cui foto sulla rivista Playboy del Novembre 1972 fu scelta (opportunamente censurata) proprio per testare algoritmi di compressione. Maggiori informazioni, oltre all’immagine stessa, le troviamo su http://www.cs.cmu.edu/ ~chuck/lennapg/. Si utilizza come   array di scalature 3 5 10 20 .

Figura 11: compressione JPEG a tratti.

37

Terzo approccio: una compressione DST. Ultimo proposito della trattazione `e quello di sostituire la DCT con un’altra trasformata discreta tra quelle presenti in letteratura: vogliamo sfruttare la DST, ossia la trasformata discreta del seno. Trattasi di un metodo di compressione puramente a fini sperimentali: vogliamo essenzialmente capire se la DST `e adatta quanto la DCT al nostro scopo o meno, cercando anche di fornire motivazioni a tale proposito. Sia v ∈ Rn , considero {s1 , . . . , sn } come base del seno:   πi sin  n + 1  2  ..   si = .   n+1 nπi  sin n+1 Le trasformate, rispettivamente diretta e inversa, sono applicazioni lineari date dalle seguenti relazioni, ricavabili analogamente a quanto fatto per la DCT. wk = (DSTn (v))k =

n X

vi sin

i=1

πki n+1

n

vk = (IDSTn (w))k =

2 X πki wi sin n + 1 i=1 n+1

Considero ora A ∈ Rm×n e voglio definire una DST matriciale. Sullo spazio delle matrici m × n, la base del seno `e {S11 , . . . , Sn1 , S21 , . . . , Smn }. ! πkj πhi 4 sin sin Sij = (n + 1)(m + 1) n+1 n + 1 h = 1, . . . , m k = 1, . . . , n Siamo pronti ad estendere le relazioni della trasformata e dell’inversa scritte in precedenza: Bhk = (DSTm×n (A))hk =

m X n X

Aij sin

i=1 j=1 m

Ahk = (IDSTm×n (B))hk =

πkj πhi sin m+1 n+1

n

XX 4 πhi πkj Bij sin sin (m + 1)(n + 1) i=1 j=1 m+1 n+1

Matlab offre una function dst vettoriale, ma non una dst2 matriciale: `e facile sfruttare la linearit` a della trasformata per fabbricarci una nostra dst2. Riducendosi al caso a noi pi` u congeniale, A ∈ Rn×n , si ha infatti: DSTn (DSTn (A)T )T = DSTn×n (A) IDSTn (IDSTn (A)T )T = IDSTn×n (A) Tutto ci` o `e svolto da queste due function.

38

% Programma che implementa una trasformata discreta del seno in % due variabili, a partire dalla function dst gi` a implementata in % Matlab. function Z=dst2(X) Y = dst(X); Z = dst(Y’)’; end % Programma che implementa una trasformata discreta inversa del seno in % due variabili, a partire dalla function idst gi` a implementata in % Matlab. function Z=idst2(X) Y = idst(X); Z = idst(Y’)’; end Osserviamo alcune peculiarit` a della trasformata del seno: anzitutto, come vediamo in Figura 12, i coefficienti elevati non andranno pi` u a posizionarsi in alto a sinistra, ma si alterneranno a coefficienti piccoli lungo tutta la superficie, come una sorta di scacchiera. Tutto questo sta a significare che dovremo rinunciare al processo di quantizzazione, o almeno non sar`a possibile utilizzare le stesse dello standard JPEG: poich`e fabbricare tali matrici `e un procedimento empirico tutt’altro che banale eviteremo categoricamente questa fase.

Figura 12: imagesc della dst2 di un blocco 8x8: verso il rosso i valori pi` u elevati. Altra caratteristica che possiamo cogliere dalla figura `e l’elevato numero di coefficienti elevati, in confronto anche alla DCT: ci`o suggerisce minore capacit`a di compressione nell’uso della DST, ancor prima di osservare qualsiasi risultato.

39

La trasformata del seno riesce a catturare con pi` u difficolt`a quelle componenti essenziali dell’immagine che la DCT raccoglie in solo 3/4 coefficienti: la severit`a di taglio non potr` a essere eccessiva, perderemmo infatti qualit`a in relativa fretta. Passiamo ora al codice che a grandi linee non differir`a molto da quello del JPEG, fatta eccezione per la fase di taglio dei valori: qua fisseremo una soglia e andremo ad eliminare quei coefficenti della DST il cui valore assoluto sar`a al di sotto della suddetta quantit`a. Per ottimizzare la computazione, su ogni blocco andremo a far uso di una matrice logica costituita da 0 e 1, con gli 1 sulle posizioni che riterremo necessarie da conservare; baster`a eseguire un prodotto termine a termine tra la DST del blocco e tale matrice e avremmo effettuato la compressione dovuta. Sceglieremo uno schema di colore in lumiocrominanza, potremo cos`ı ricondurci con facilit` a a quanto fatto in precedenza e allo stesso tempo avremo un confronto pi` u diretto. Utilizziamo i medesimi programmi in C crea matrici.c e crea ppm.c, oltre alle function di conversione dall’RGB all’YCbCr e viceversa. Riportiamo di seguito la function di compressione, quella di decompressione e lo script .sh di esecuzione. % Programma che esegue una compressione mediante DST % sulle componenti YCbCr, eliminando i valori aventi % modulo sotto la soglia data. % k rappresenta la dimensione di ciacun sottoblocco in % cui viene divisa la foto. % z rappresenta il fattore di normalizzazione. function compress_dst(k,soglia,z) % Prima parte: adattamento immagine alla suddivisione in blocchi. disp(’Adatto componenti...’) dim = load(’dimensione_foto’,’-ascii’); largh = dim(1); alt = dim(2); foto_Y = load(’matr_Y’,’-ascii’); foto_Cb = load(’matr_Cb’,’-ascii’); foto_Cr = load(’matr_Cr’,’-ascii’); foto_Y = reshape(foto_Y,largh,alt)’; foto_Cb = reshape(foto_Cb,largh,alt)’; foto_Cr = reshape(foto_Cr,largh,alt)’; if (mod(largh,k) ~= 0) % Adattamento larghezza foto. for i=1:k-mod(largh,k) foto_Y(1:end,largh+i) = 128; foto_Cb(1:end,largh+i) = 128; foto_Cr(1:end,largh+i) = 128; end end if (mod(alt,k) ~= 0) % Adattamento altezza foto. 40

for i=1:k-mod(alt,k) foto_Y(alt+i,1:end) = 128; foto_Cb(alt+i,1:end) = 128; foto_Cr(alt+i,1:end) = 128; end end foto_Y = foto_Y(1:alt+k*(mod(alt,k) ~= 0)-mod(alt,k),1:largh+ +k*(mod(largh,k) ~= 0)-mod(largh,k)); foto_Cb = foto_Cb(1:alt+k*(mod(alt,k) ~= 0)-mod(alt,k),1:largh+ +k*(mod(largh,k) ~= 0)-mod(largh,k)); foto_Cr = foto_Cr(1:alt+k*(mod(alt,k) ~= 0)-mod(alt,k),1:largh+ +k*(mod(largh,k) ~= 0)-mod(largh,k)); sz = size(foto_Y); dim(1) = sz(2); dim(2) = sz(1); % Seconda parte: calcolo della DST e compressione. disp(’Calcolo DST e comprimo...’); matr_Y = zeros((dim(1)/k)*(dim(2)/k),k); matr_Cb = zeros((dim(1)/k)*(dim(2)/k),k); matr_Cr = zeros((dim(1)/k)*(dim(2)/k),k); cursore_matr = 1; for i=1:k:dim(1)-k+1 for j=1:k:dim(2)-k+1 trasf_Y = dst2(foto_Y(j:j+k-1,i:i+k-1)); % Calcolo trasformate. trasf_Cb = dst2(foto_Cb(j:j+k-1,i:i+k-1)); trasf_Cr = dst2(foto_Cr(j:j+k-1,i:i+k-1)); trasf_Y = trasf_Y .* (abs(trasf_Y) > soglia); % Elimino valori sotto la soglia data. trasf_Cb = trasf_Cb .* (abs(trasf_Cb) > soglia); trasf_Cr = trasf_Cr .* (abs(trasf_Cr) > soglia); matr_Y(cursore_matr:cursore_matr+k-1,:) = floor(trasf_Y*z); matr_Cb(cursore_matr:cursore_matr+k-1,:) = floor(trasf_Cb*z); matr_Cr(cursore_matr:cursore_matr+k-1,:) = floor(trasf_Cr*z); cursore_matr = cursore_matr+k; end end disp(’Scrivo i dati...’); matr_Y = sparse(matr_Y); % Acquisizione come matrice sparsa. matr_Cb = sparse(matr_Cb); matr_Cr = sparse(matr_Cr); save(’compressa.mat’,’matr_Y’); save(’compressa.mat’,’matr_Cb’,’-append’); save(’compressa.mat’,’matr_Cr’,’-append’); end

41

# Script di esecuzione: alg. di compressione DST. #!/bin/bash echo "Soglia di compressione?" read SOGLIA echo "Creo matrice per Octave/Matlab..." gcc -o crea_matrici crea_matrici.c ./crea_matrici echo "RGB_to_YCbCr();compress_dst(8,${SOGLIA},20)" | matlab echo "Elimino file non piu’ utili..." rm "matr_R" rm "matr_G" rm "matr_B" rm "matr_Y" rm "matr_Cb" rm "matr_Cr" echo "Fatto!" % Programma che esegue la decompressione calcolando per ogni blocco la trasformata % del coseno inversa a partire dalla matrice sparsa dei coefficienti conservati. % k rappresenta la suddivisione dei sottoblocchi (k=8 per il JPEG). % z rapprensenta il fattore di normalizzazione. function decompress_dst(k,z) load(’compressa.mat’); %Carico elementi necessari. matr_Y = full(matr_Y); %Recupero della forma completa. matr_Cb = full(matr_Cb); matr_Cr = full(matr_Cr); dim_foto = load(’dimensione_foto’,’-ascii’); largh = dim_foto(1)+k*(mod(dim_foto(1),k) ~= 0)-mod(dim_foto(1),k); alt = dim_foto(2)+k*(mod(dim_foto(2),k) ~= 0)-mod(dim_foto(2),k); cursore_matr = 1; foto_Y = zeros(alt,largh); %Inizializzo matrice dell’immagine. foto_Cb = zeros(alt,largh); foto_Cr = zeros(alt,largh); blocchi_largh = largh/k; blocchi_alt = alt/k; disp(’Decomprimo la foto...’); for i=1:blocchi_largh for j=1:blocchi_alt %Calcolo l’antitrasformata del seno. antitrasf_Y = round(idst2(matr_Y(cursore_matr:cursore_matr+k-1,:)/z)); antitrasf_Cb = round(idst2(matr_Cb(cursore_matr:cursore_matr+k-1,:)/z)); antitrasf_Cr = round(idst2(matr_Cr(cursore_matr:cursore_matr+k-1,:)/z)); foto_Y((j-1)*k+1:j*k,(i-1)*k+1:i*k) = antitrasf_Y; foto_Cb((j-1)*k+1:j*k,(i-1)*k+1:i*k) = antitrasf_Cb; 42

foto_Cr((j-1)*k+1:j*k,(i-1)*k+1:i*k) = antitrasf_Cr; cursore_matr = cursore_matr+k; end end foto_Y = foto_Y(1:dim_foto(2),1:dim_foto(1)); foto_Cb = foto_Cb(1:dim_foto(2),1:dim_foto(1)); foto_Cr = foto_Cr(1:dim_foto(2),1:dim_foto(1)); foto_Y = foto_Y’; %Calcolo trasposta, in modo da scrivere per riga. foto_Cb = foto_Cb’; foto_Cr = foto_Cr’; Y = fopen(’matr_Y’,’w’); %Salvo matrice. Cb = fopen(’matr_Cb’,’w’); Cr = fopen(’matr_Cr’,’w’); fprintf(Y,’%d\n’,foto_Y); fprintf(Cb,’%d\n’,foto_Cb); fprintf(Cr,’%d\n’,foto_Cr); fclose(Y); fclose(Cb); fclose(Cr); end

# Script di esecuzione: alg. di decompressione DST. #!/bin/bash echo "decompress_dst(8,20);YCbCr_to_RGB()" | matlab echo "Creo file .pgm..." gcc -o crea_ppm crea_ppm.c ./crea_ppm echo "Elimino files non piu’ utili..." rm "matr_Y" rm "matr_Cb" rm "matr_Cr" rm "matr_R" rm "matr_G" rm "matr_B" echo "Fatto!" Effettuiamo ora qualche prova sul programma: considerando ancora una volta l’immagine della Canon (49,7 MB) facciamo variare le soglie di compressione al fine di testare l’effettiva capacita dell’algoritmo DST. Fissando il valore di soglia a 20 non si notano imperfezioni e si riesce a ridurre lo spazio a 7,9 MB; con soglia di taglio 60 si scende a 4,3 MB, con 100 a 3,2 MB, con 200 a 1,9 MB. Da notare che queste ultime tre prove riportano difetti visivi rispetto all’immagine originale: la carenza di qualit`a si inizia a notare dal perimetro di ogni suddivisione e va ad allargarsi man mano che il taglio `e pi` u severo, si ha una visione delle suddivisioni 8x8 sempre pi` u netta.

43

Figura 13: l’immagine intera non compressa (.ppm)

Figura 14: particolari con soglie di taglio 20, 60, 100, 200.

44

Conclusioni Siamo giunti alla fine del percorso e vediamo di tirare il punto della situazione: abbiamo effettivamente confrontato un algortimo commerciale quale il JPEG con due tentativi pi` u improvvisati. Se volessimo stilare una sorta di classifica non avremmo ovviamente problemi su chi posizionare al primo posto: un poco pi` u difficile, ma neppur troppo, stabilire invece il migliore tra l’approccio ` ovvia la superiorit`a della trasformata del seno in termini di SVD e il DST. E compressione: da un lato abbiamo conservato solo alcuni coefficienti numerici e dall’altro diversi vettori singolari, una mole di dati maggiore. In termini di qualit` a invece l’SVD non ha nulla a che invidiare agli altri; grazie all’apporto dato dalla suddivisione dell’immagine questo algoritmo lavora egregiamente sulle variazioni non drastiche, mentre presenta dei problemi negli alti cambi di tono solo in caso di una conservazione di valori singolari troppo esigua. Ultima puntualizzazione che mi pare necessaria: il codice che abbiamo scritto non rappresenta il massimo dell’efficienza in termini di tempo di esecuzione, ci proponiamo di ”volare basso” e considerare tutto ci`o mera presentazione di esempi e di principi, di confronti tra possbili idee matematiche atte a suggerire metodi di compressione di immagini digitali.

45