G. Naldi, T. Scapolla

Sunto. In questo lavoro si vogliono descrivere alcuni aspetti relativi all’utilizzo di recenti tecniche per la compressione e l’analisi di immagini nell’ambito della procedura investigativa che utilizzi le

impronte digitali. I procedimenti considerati per la trasformazione e/o compressione di una immagine si basano sulla trasformata wavelet e sulla compressione frattale. Si considererà inoltre l’accoppiamento tra sistema visivo umano e compressione wavelet e si descriveranno possibili sviluppi futuri per procedimenti automatizzati di identificazione di impronte digitali .

Le creste cutanee sono rilievi determinati in alcune regioni del corpo (palma delle mani, pianta dei piedi, polpastrelli delle dita) dalla sporgenza delle papille dermiche allineate (scoperte da M. Malpighi nel 1664). Le creste cutanee hanno la larghezza di 0.2-0.7 mm. e sono visibili a occhio nudo. Il complesso delle varie creste cutanee forma dei caratteristici disegni studiati fin dal 1823 dal fisiologo J.H. Purkinje e in seguito da F. Galton (Finger prints , Londra 1892). Essendo tali disegni in rilievo, possono lasciare impronte sulle superfici liscie. Dato che i disegni digitali e palmari sono strettamente individuali, le relative impronte possono essere utilizzate nel processo di identificazione nei servizi di polizia. La dattiloscopia è lo studio scientifico e pratico di questi disegni digitali e delle relative impronte. Le impronte digitali (cioè del polpastrello delle dita delle mani) presentano differenze generali (utili per una loro classificazione) e differenze particolari (utili per l’identificazione individuale). Queste caratteristiche stanno alla base delle due principali applicazioni pratiche che sono state fatte : l’identificazione dei recidivi, l’individuazione delle impronte lasciate sul luogo del reato.

Per l’identificazione si raccolgono le impronte delle dieci dita sulla scheda segnaletica di ogni individuo e poi si ripartiscono le schede in appositi casellari. La distribuzione si fa in base alla formula dattiloscopica che varia da sistema a sistema. Il sistema italiano, dovuto al Gasti con succesive integrazioni del Sorrentino, prevede più categorie. L’identità si dimostra col constatare che due impronte appartengonom allo stesso tipo, con lo stesso numero di linee, col rilievo di un certo numero di punti caratteristici (12-15 punti caratteristici) egualmente situati (biforcazioni, interruzioni, anelli, ecc.) e con l’assenza di punti caratteristici in contrasto.

Per velocizzare ed agevolare le operazioni di identificazione e di classificazione delle impronte digitali, diventa necessario prevedere la memorizzazione delle impronte stesse su supporto elettronico tramite un opportuno formato ‘‘digitalizzato’’ . Per esempio, come fatto in America dall’FBI, ogni impronta su scheda, potrebbe essere digitalizzata in una immagine i bianco e nero con scala di grigio a 8 -bit e con una risoluzione di 500 dpi (dots per inch). In questo modo ogni scheda richiederebbe circa 10 megabytes per una opportuna memorizzazione. Data la quantità di impronte (schede) da memorizzare ed il numero di nuove acquisizioni giornaliere ( sempre con riferimento all’FBI, si hanno circa trentamila nuove acquisizioni per giorno, ovverosia 300 gigabytes) diventa essenziale un processo di compressione di queste immagini.

In Figura 1 si mostra un esempio di impronta digitale ‘‘digitalizzata’’ per la memorizzazione su un opportuno supporto elettronico. Questo esempio è stato tratto dalla biblioteca di immagini disponibili con il Toolbox Wavelet di Matlab.

Una metodologia generale per la compressione delle immagini consiste nel considerare una opportuna rappresentazione del segnale (preferibilmente tramite una ‘‘trasformata rapida’’) che permetta di utilizzare un numero minore di bit per la memorizzazione mantenendo una accuratezza accettabile (eliminando le informazioni ‘‘rindondanti’’). Uno degli approcci più utilizzati impiega l’analisi di Fourier e le relative trasformate rapide per lo sviluppo del segnale in basi trigonometriche (funzioni seno e coseno a frequenze differenti). In particolare, un diffuso metodo di codifica è lo standard JPEG (Joint Photographic Expert Group).

Il processo di trasformazione di JPEG prevede diversi passi. Inizialmente l’immagine viene suddivisa in blocchi contigui di 8´ 8 pixels. Ogni blocco è quindi rappresentato da una matrice di 8 righe e 8 colonne di numeri interi variabili nell’intervallo [0, 255]. Quindi si normalizzano tali valori in modo tale che varino tra -128 e 127. Nel secondo passo si applica la trasformata coseno discreta ad ogni blocco ottenendo una nuova matrice trasformata, sempre di dimensioni 8´ 8. La matrice trasformata, che contiene le componenti per differenti frequenze, viene ‘‘quantizzata’’ suddividendo le componenti, opportunamente normalizzate, in sottointervalli e stabilendo un valore intero rappresentativo per ogni sottointervallo. Questo passo è essenziale per stabilire il livello di compressione e la qualità dell’immagine compressa. La matrice quantizzata viene memorizzata in un vettore tramite una lettura a ‘‘zig-zag’’ della matrice stessa. Infine il vettore viene codificato tramite una codifica run-length e, successivamente, tramite un metodo entropico quale il metodo di Huffman. Il procedimento di decompressione è facilmente ottenibile ‘‘invertendo’’ i singoli passi utilizzati nella fase di compressione.

Generalmente il rapporto di compressione CR è definito da

CR = __________________________________________

Utilizzando JPEG, anche per moderati valori del livello di compressione, ci si accorge che si formano degli artefatti nell’immagine decompressa dovuti all’interfacciamento dei singoli blocchi di 8 pixels per lato. Nel caso delle impronte digitali questo fenomeno non è accettabile dato che può influenzare il processo di identificazione. Una discussione di questo aspetto si può trovare in C.M. Brislawn , mentre ulteriori informazioni si possono ricavare su un link World Wide Web appositamente creato : http ://www.c3.lanl.gov/~ brislawn/main.html. In figura 2 si mostra un tipico risultato (qui preso da C.M. Brislawn) ottenibile con lo standard JPEG ed una memorizzazione che utilizzi 0.6 bits/pixel.

Recentemente le basi wavelets hanno fornito una importante alternativa sia per la compressione delle immagini che nei più svariati ambiti applicativi : dall’edge detection agli operatori integrali, dall’analisi di strutture frattali agli algoritmi adattativi nella discretizzazione di equazioni differenziali, dall’analisi di segnali biomedici alla simulazione di strutture molecolari.

In particolare tali basi permettono di effettuare una analisi di Fourier localizzata permettendo l’individuazione delle componenti ‘‘importanti’’ di un segnale (funzione, immagine) sia nello spazio che in frequenza.

Il termine wavelet (originalmente chiamata ‘‘wavelet of constant shape’’) è stato introdotto da J. Morlet ed in generale denota una funzione y(x) definita sulla retta reale R che, soggetta alle operazioni di traslazione di un intero k e dilatazione di un fattore 2^j ( j intero), fornisce una base ortonormale per lo spazio L²(R) delle funzioni a quadrato integrabile. In altre parole l’insieme di funzioni :

 { y_jk (x) = 2^j/2 y (2^j x -k), j,k Î Z } (1)

forma un sistema ortonormale completo per L²(R ).

In figura 3 si mostrano i grafici di alcune funzioni wavelet, la wavelet di Haar (introdotta da A. Haar nel 1910 in altro ambito) può essere storicamente considerata la prima wavelet.

La base (1) permette di ottenere una decomposizione gerarchica di un dato segnale, questa proprietà è stata formalizzata nel concetto di Multiresolution Analysis da S. Mallat e Y. Meyer. La decomposizione in una base di wavelets di un segnale può essere efficientemente fatta tramite una trasformata rapida (con una complessità linearmente proporzionale alla lunghezza del segnale discreto). Inoltre la decomposizione che si ottiene permette di individuare solo le componenti essenziali del segnale in esame grazie alle proprietà di localizzazione spazio/frequenza delle funzioni wavelets.

Un altra tecnica introdotta recentemente per la compressione delle immagini utilizza idee provenienti da una teoria matematica nota come iterated function systems (IFS). Tale tecnica sfrutta l’autosimilità di parti dell’immagine e la potenzialità dei frattali per la modellizazione dei fenomeni naturali. La teoria è stata inizialmente sviluppata da J. Hutchinson nel 1981, in seguito M. Barnsley ed il suo gruppo di ricerca presso il Georgia Institute of Technology, hanno incominciato a sfruttare l’approccio IFS per modellizare strutture naturali quali nuvole, alberi, foglie. Nel 1989 A. Jacquin, un allievo di Barnsley, realizzò il primo sistema automatico di codifica frattale di una immagine.

Il codice forniva buoni indici di compressione ma risultava ‘‘lento’’ ; miglioramenti e varianti sono state proposte in seguito. In particolare è risultato fondamentale il lavoro di Y. Fisher che ha utilizzato uno schema di codifica frattale adattivo con una struttura dati tipo quadtree per i blocchi che compongono l’immagine da comprimere.

In questo ultimo approccio infatti si partiziona l’immagine in vari blocchi (che si possono anche sovrapporre) e per ogni blocco si cerca un opportuno IFS che lo possa rappresentare con una adeguata accuratezza. In altri termini si stabilisce una porzione di figura (range block R) ed una parte più grande (domain block D) ed una trasformazione affine T tale che T(D), o una sua iterata Tⁿ, sia una buona approssimazione del blocco R. Solo la trasformazione sarà codificata contenendo tutte le informazioni necessarie per ricostruire la porzione di immagine selezionata. Naturalmente l’applicazione T deve soddisfare opportune proprietà, in particolare deve essere una applicazione contrattiva (in caso contrario la decodifica non potrà essere fatta).

I risultati ottenuti con la tecnica di compressione frattale sembrano promettenti anche se il processo di codifica risulta essere computazionalmente più oneroso mentre quello di decodifica più semplice.

Consideriamo ora la nozione fondamentale di Analisi Multirisoluzione (nel seguito abbreviata con MRA). Rimandiamo naturalmente ai riferimenti bibliografici indicati alla fine del presente lavoro per una trattazione più approfondita della teoria delle wavelets.

Consideriamo il problema della costruzione di un sistema ortonormale completo per lo spazio L²(R ) partendo dalla funzione caratteristica per l’intervallo unitario [0,1) :

 c _[0,1) (x) = f (x) (2)

E’ chiaro che l’insieme { f (x - k ), k Î Z } è un sistema ortogonale in L²(R ) ma non è completo dato che genera il solo sottospazio chiuso V ₀ delle funzioni costanti a tratti con ‘‘salti’’ possibili solo per valori interi, per esempio la funzione c _[0,1/2) (x) non appartiene a tale sottospazio.

Per includere più funzioni consideriamo allora le versioni dilatate, con fattore di dilatazione 2^j e j intero, traslate di un valore k intero e normalizzate, della funzione caratteristica f (x) e:

 f _jk (x) = 2^j/2 f (2^j x - k) (3)

Per un valore fissato del parametro intero j l’insieme delle funzioni (3) genera il sottospazio lineare chiuso denotato con V_j . Dato che ogni funzione di fÎ L²(R ) può essere ‘‘arbitrariamente’’ approssimata da una funzione costante a tratti f_j con salti nei razionali diadici del tipo x_jk=k/2^j , ne segue che l’unione È _j V_j è densa in L²(R ). Da quanto detto si deduce che il sistema delle funzioni (3) al variare degli indici interi j,kÎ Z è un sistema completo ma si è persa l’ortogonalità (per esempio le funzioni caratteristiche rappresentate da f (x) e f (2x) non sono ortogonali). Per ovviare a ciò si può introdurre una nuova funzione y (x) appartenente a V₁ ed ortogonale alla funzione f (x). Un semplice calcolo mostra che la funzione y (x)= f (2x)- f (2x-1) soddisfa ad entrambe le richieste. Si deduce quindi che l’insieme delle funzioni y _jk (x) = 2^j/2 y (2^j x - k) al variare di j e k negli interi formano un sistema ortonormale completo : il sistema di Haar. La funzione f è nota come funzione scaling di Haar mentre la funzione y mother wavelet di Haar (si veda la figura 3). Nel sistema di Haar sono validi i seguenti sviluppi in serie, sempre nell’ambito dello spazio L²(R ) (qui < × ,× > indica il prodotto scalare in tale spazio, <f,g >=ò R f(x) g(x)dx) :

 f(x) = S _{jS k} d_jk y _jk (x) , d_jk =<f, y _jk > (4.a)

f(x) = S _k c_mk f _mk + S _{j³ m S k} d_jk y _jk (x) , d_jk =<f, y _jk > , c_mk =<f, f _mk > (4.b)

a cui corrispondono diverse decomposizioni in sottospazi. Formalizziamo tale costruzione ad un caso più generale.

Definiamo MRA ortogonale di L²(R ) una successione { V_j} di sottospazi chiusi a cui sia associata una funzione f Î V₀ (detta funzione scaling) e tali che

(i) l’insieme delle funzioni { f (x-k), kÎ Z } è una base ortonormale per V₀ ;

(ii) i sottospazi formano una successione crescente : .... Ì V_-1 Ì V₀ Ì V₁ Ì ... ;

(iii) cambiamento di scala f(× ) Î V_j Û f(2× ) Î V_j+1 ; invarianza per traslazioni f(× )Î V₀ Û f(× -k)Î V₀

Chiaramente l’insieme delle funzioni (3) per j fissato e k variabile negli interi formano una base ortonormale per lo spazio V_j . Usualmente si fanno richieste di regolarità e di decadimento (comportamento all’infinito della funzione e di alcune derivate successive) della funzione f .

Dalla funzione scaling f si costruisce poi la funzione wavelet y (solitamente tale costruzione utilizza la trasformata di Fourier ) cioè una funzione y tale che l’insieme { y (x-k)} formi una base ortonormale del sotospazio W₀ complemento ortogonale di V₀ in V₁, V₁=V_0Å W₀. Dalla definizione della MRA si deduce che dal sottospazio W₀ possiamo costruire valgono le seguenti decomposizioni dello spazio delle funzioni a quadrato integrabile:

L²(R )=Å j W_j , L²(R )=V_m Å _{j³ m} W_j

Si osserva che esistono delle relazioni funzionali ricorsive sia per la funzione f che per la funzione y , infatti, dato che f , y Î V₁ devono esistere dei coefficienti { C_k} e { D_k} tali che :

f (x)= S _k Ö 2 C_k f (2x-k), (5.a)

y (x)= S _k Ö 2 C_k f (2x-k) (5.b)

Inoltre se, almeno localmente, i polinomi fino ad un certo grado appartengono al sottospazio V₀ (questo proprietà influenza la qualità dell’approssimazione che si ottiene sostituendo ad una data funzione f la proiezione f_j su un opportuno sottospazio V_j), la funzione y avrà un certo numero di momenti nulli :

ò R y (x) x^k =0 , k=0,1,...,r. (6)

Consideriamo una funzione fÎ V₁ possiamo rappresentarla come :

f(x)=S _k a_1k f _1k = _{S k} a_{1k Ö} 2 _f (2x-k) (7)

ma dalla decomposizione di V₁ nei sottospazi V₀ e W₀ si ottiene anche

f(x)= S _k a_0k f (x-k) + _{S k} b_0k y (x-k) (8)

Sostituendo le (5.a) e (5.b) nella (8) ed eguagliando i coefficienti corrispondenti alla stessa funzione f _1k (abbiamo basi ortonormali) si ottiene :

a_1m=S _k C_m-2k a_0k + _{S k} D_m-2k b_0k (9)

 che rappresenta il passo di ricostruzione da V₀, W₀ a V₁. Per la parte di decomposizione basta osservare che, tenendo conto sempre di (5.a) e (5.b) :

b_0m=ò R y (x-k) f(x) dx = S _k D_k-2m a_1k (10.a)

a_0m=ò R f (x-k) f(x) dx = S _k C_k-2m a_1k (10.a)

Nelle sommatorie in (9), (10.a), (10.b) si riconosce la struttura algebrica dei prodotti di convoluzione ed in effetti corrispondono alla applicazione di opportuni filtri discreti alle successione dei coefficienti { a_jk } e { b_jk } (immaginando che i procedimenti di decomposizione/ricostruzione siano iterati).

Queste osservazioni hanno portato ad un collegamento naturale (dal punto di vista dell’implementazione) tra decomposizione gerarchica in una MRA ed i banchi di filtri utilizzati in ingegneria, inoltre pongono le basi per lo sviluppo di algoritmi rapidi.

In figura 4 si mostra una decomposizione su più scale di una data funzione, grafico con etichetta ‘originale’, cioè delle componenti nei sottospazi V_j e W_j. In questo esempio si è utilizzata la MRA di Haar, quindi i coefficienti C e D sono rispettivamente { 1/Ö 2, 1/Ö 2} , { 1/Ö 2, -1/Ö 2} .

Si può notare comunque euristicamente che le componenti negli spazi W_j sono piccole dove la funzione è regolare ed è quindi approssimata con sufficiente accuratezza dalla proiezione su V_j. Nel caso considerato questa ultima proiezione è facilmente ottenibile tramite medie integrali locali della funzione per suddivisioni con ampiezza decrescente. In effetti le componenti nei W_j sono dette dettagli e codificano l’informazione necessaria per passare dalla descrizione ad una certa scala della funzione (segnale, immagine,...) alla descrizione alla scala successiva. In tali dettagli si codificano quindi eventuali irregolarità locali della funzione che si sta analizzando.

Esistono diverse famiglie di wavelets (a supporto compatto, interpolanti, biortogonali,...) e, per ogni famiglia, si hanno differenti membri con caratteristiche analoghe ma differenti parametri nei coefficienti C e D. La trasformata discreta, relazioni (9), (10.a), (10.b), ha in ogni caso la medesima struttura. La scelta di una particolare famiglia per un determinato problema non sempre è semplice, in ogni caso nell’ambito della compressione delle immagini la sperimentazione effettuata indica nelle wavelet biortogonali una scelta ottimale.

Come sono modellate matematicamente le immagini ? Solitamente un immagine formata da una griglia di pixels, ognuno con un livello di grigio, è vista come una funzione z=f(x,y) in due variabili reali (x,y)Î I´ J, dove I=[a,b] e J=[c,d] sono due intervalli limitati. Quando ci si riferisce ad una immagine ci si riferisce quindi ad una funzione che per ogni punto (x,y) fornisce il corrispondente livello di grigio. Nel caso si voglia considerare la ‘‘differenza’’ tra due immagini (per esempio tra immagine originale ed immagine compressa) occorre definire anche una metrica e quindi una funzione che permetta di stabilire la ‘‘distanza’’ tra due immagini. Due esempi di funzioni distanza tra l’immagine f e l’immagine g sono dati dalle seguenti definizioni :

d8 (f,g)=sup_{(x,y)Î I´ J} | f(x,y) -g(x,y) | , d₂(f,g)=( ò _{I´ J} (f (x,y) - g(x,y) )² dx dy )½

Negli esperimenti del seguito utilizzeremo la versione discreta della distanza d₂ anche se il problema della metrica ‘‘giusta’’ nella teoria matematica delle immagini resta un problema aperto. In ogni caso la scelta può essere influenzata dalla particolare applicazione che si sta considerando. Per esempio nell’ambito della segmentazione delle immagini e dell’edge detection una buona scelta sembra essere la variazione totale dell’immagine, ò |Ñ f(x,y)|.

Per effettuare la trasformata wavelet di una immagine discreta I (quindi partendo da una matrice di valori) si possono considerare le trasformate unidimensionali delle righe ottenendo la versione I_r , quindi si effettuano le trasformate unidimensionali delle colonne di I_r ottenendo l’immagine trasformata I_W. I dettagli saranno quindi organizzati in sottomatrici di I_W .

Per utilizzare questa struttura gerarchica ai fini della compressione delle immagini (per altre applicazioni si veda la bibliografia al termine di questo articolo) si selezionano proprio i dettagli.

Naturalmente si perderanno alcune informazioni relative al segnale ma, nelle applicazioni, si sfrutta la naturale correlazione tra i valori dei pixels dell’immagine stessa.

Per costruire una MRA in L²(R ²) e le relative funzioni base si effettua il prodotto tensoriale delle basi unidimensionali in x ed y. Abbiamo quindi una funzione scaling, che chiameremo ancora f (x,y)= f (x) f (y) e tre funzioni wavelets y ¹= _f (x) y (y), y ²= _f (y) y (x), y ³ = _y (x) y (y), in un certo senso corrispondenti ai dettagli orizzontali, verticali e ,rispettivamente, obliqui. Dal punto di vista discreto ogni sottomatrice corrispondente ad un sottospazio V_j bidimensionale verrà suddivisa, nella decomposizione, in quattro sottomatrici corrispondenti alle quattro funzioni f , y ¹, y ², y ³ _.

Indichiamo con f_I la funzione corrispondente all’immagine I, per funzioni in due valgono le decomposizioni analoghe alle (4.a) e (4.b), anche se ora abbiamo in un certo senso più funzioni wavelets :

f(x,y)= S _k a_0k f (x,y) +S _p S _{j S k} b_pjk y ^p_jk(x,y), (11)

dove p=1,2,3. Solitamente i coefficienti a_0k vengono mantenuti (nella terminologia dei filtri rappresenta la banda base a cui aggiungeremo i dettagli succesivi) mentre possiamo selezionare i coefficienti b_pjk . Un aspetto importante deve essere considerato per l’effettiva implementazione della trasformata wavelet : il trattamento del bordo dell’immagine. Infatti fino ad ora abbiamo trattato il caso di funzioni definite su tutta la retta reale o su tutto il piano R ². Se la funzione, come accade nelle applicazioni, è definita solo su un particolare sottoinsieme, per esempio, del piano occorre utilizzare una opportuna estensione dell’immagine stessa vicino ai bordi.

Una strategia adatta alla nostra applicazione consiste nell’estensione simmetrica del segnale in modo tale da non introdurre discontinuità artificiali al bordo, discontinuità che influenzerebbero i coefficienti relativi alle funzioni wavelets in (11).

Indichiamo con S l’insieme delle triple (p,j,k) degli indici selezionati secondo una determinata strategia. Sono possibili varie richieste, per esempio :

-- fissare la cardinalità N dell’insieme S (massima memoria utilizzabile) ;

-- stabilita una norma ║ × ║ cercare S che permetta di minimizzare l’errore ║ f-g║ , dove g=S _k a_0k f (x,y) +S _p S _{j S k} b_pjk y ^p_jk(x,y), (j,k,p)Î S, oppure fissato il valore dell’errore determinare un possibile insieme S ;

-- fissare il valore minimo di | b_pjk| con (j,k,p)Î S.

Nell’applicazione che stiamo analizzando si effettua una selezione dei coefficienti b_pjk tramite opportuni valori di soglia l (j,e , ║ × ║ ) che possono dipendere, a priori, dal livello di decomposizione j, da una prefissata tolleranza e dalla norma utilizzata per la valutazione dell’errore. La scelta più semplice consiste nell’ hard thresholding : l è costante e se | b_pjk | ³ l il coefficiente verrà mantenuto altrimenti sarà posto uguale a zero (quindi eliminato).

Prima di tale selezione si effettua una normalizzazione dei coefficienti b_pjk in funzione della scala j e del valore max_k |b_pjk| per ogni valore dell’indice p.

Riportiamo nella tabella 1 che segue la percentuale di elementi ’eliminati’ con tale processo di selezione utilizzando l’impronta digitale di figura 1 ed wavelet biortogonali uniti alla normalizzazione dei coefficienti prima del confronto con la soglia l .

In figura 5 invece mostriamo l’impronta digitale ottenuta tramite lo sviluppo approssimato dopo la selezione dei coefficienti con l ==0.1.

Osserviamo che nel caso delle wavelet biortogonali si usano differenti funzioni per la fase di decomposizione e per la fase di ricostruzione. In questo modo si possono richiedere più proprietà alla funzione wavelet y che invece risulterebbero incompatibili nel caso semplicemente ortogonale : supporto compatto (filtri di lunghezza finita), simmetria (fase lineare dei filtri), filtri discreti a valori nei razionali (rapidità di calcolo).

Nel caso di hard thresholding si introduce una discontinuità nella funzione di selezione per i valori ± l , per ovviare a ciò si può utilizzare il soft thresholding in cui, se indichiamo con b^S_pjk il valore dei coefficienti dettaglio dopo il processo di selezione, abbiamo b^S_pjk=0 se | b^S_pjk|£ l e b^S_pjk=sign(b^S_pjk)(| b^S_pjk|-l ) altrimenti. I metodi di thresholding possono anche essere utilizzati per problemi di denoising dell’immagine fornendo quindi contemporaneamente un modo per decorrelare le informazioni e per togliere eventuale rumore dall’immagine stessa.

Un passo importante consiste ora nella quantizzazione dei coefficienti ottenuti dopo l’operazione di thresholding, chiamiamo tali coefficienti b^T_pjk. In generale tali coefficienti sono reali (floating point) e si vuole sostituirli con un numero intero : in un certo senso il passo di quantizzazione equivale a fare una ‘‘sofisticata’’ operazione di troncamento. In figura 6 si mostra un esempio di funzione di quantizzazione in cui ad ogni intervallo [a _m,b _m) dei valori dei coefficienti b^T_pjk già selezionati si fa corrispondere un valore intero.

Per ora nell’applicazione qui considerata si è utilizzata per ogni scala una suddivisione uniforme dell’intervallo [max_kb^T_pjk , min_k b^T_pjk].

Sottolineiamo che il passo di quantizzazione è sicuramente il passo più ‘‘delicato’’ ed importante ai fini della compressione dell’immagine. Nello sviluppo futuro di considererà invece una quantizzazione che prenda in considerazione la particolare statistica per ogni scala e tipo di dettaglio.

Notiamo anche che mentre dal punto di vista logico è utile tenere separate le fasi di selezione con thresholding e di quantizzazione, dal punto di vista dell’implementazione del metodo di compressione i due processi possono essere fatti contemporaneamente tramite una opportuna variante del metodo di hard thresholding. Osserviamo anche che i coefficienti relativi alla funzione scaling f non vengono sottoposti ai procedimenti di thresholding/quantizzazione in quanto si pensa che il contenuto informativo di tali coefficienti è essenziale e non può essere ‘‘trattato’’ come dettaglio.

Il passo finale dell’algoritmo di compressione consiste in una ulteriore codifica dei coefficienti quantizzati b^Q_pjk che sono ora rappresentati tramite un vettore di interi. Tale codifica finale utilizzerà un metodo di compressione lossless, cioè senza perdita di informazioni, basato sull’algoritmo di Huffman. Globalmente il metodo di compressione fornisce fattori di compressione di 0.6-0.9 bit per pixel con una buona qualità per quanto riguarda l’immagine ricostruita. Per ottenere quest’ultima immagine si invertono i singoli passi visti sopra.

Abbiamo accennato al problema della scelta di una metrica opportuna per valutare la distanza tra due immagini (e quindi valutare anche la qualità dell’immagine compressa). Dal punto di vista discreto dobbiamo stabilire una funzione distanza tra due matrici di valori (che rappresenteranno il livello di grigio dei singoli pixel). Supponiamo che le immagini siano rappresentate da matrici di n righe ed m colonne, gli indicatori di errore più frequentemente utilizzati sono i seguenti  (f e g sono le due matrici):

-- errore quadratico medio MSE : S _i S _j ( f_ij - g_ij )² / (m´ n)

-- errore quadratico normalizzato NRMSE : 100 ´ (MSE)^1/2 / 255

-- rapporto segnale/rumore  PSNR : 10 log₁₀ ( 255²/MSE)

-- errore massimo PE : max_ij | f_ij - g_ij |

-- errore massimo normalizzato NPE : 100´ PE/255

Le stime PE ed NPE considerano differenze locali mentre l’errore quadratico medio è una stima globale ma non risulta adatto per problemi in cui eventuali caratteristiche dell’immagine quali gli edges. Nell’ambito della compressione delle immagini differenti approcci possono invece essere introdotti per valutare l’errore di ricostruzione considerando le proprietà fisiologiche dell’occhio umano. In particolare, sotto opportune condizioni sperimentali di isotropia, linearità del sistema visivo, luminosità, si è constato che l’occhioumano ‘‘pesa’’ in modo differente frequenze spaziali differenti. La forma di trasferimento h(r) è stata determinata sperimentalmente ed il grafico qualitativo è mostrato in figura 7.

Sembra quindi ragionevole, utilizzando la trasformata di Fourier, introdurre tale caratteristica nella valutazione dell’errore. Indichiamo con I(x,y) ed I^R(x,y) , rispettivamente l’immagine originale e l’immagine ricostruita, e con i(h ,V ), i^R(h ,V ) le loro trasformate di Fourier. L’errore, pesato tramite la funzione h è definito da

|| I - I^R ||²_h = ò ò R ´ R h( (h ² + V ²)^1/2) | i(h ,V ) - i^R (h ,V )|² dh dV (12)

Partendo da una MRA ortogonale o biortogonale è possibile introdurre una nuova MRA che consideri come prodotto scalare non il classico prodotto in L²(R ) ma il prodotto (12). In questa costruzione però risulta esserci una dipendenza dei filtri dal livello j di decomposizione, non essendoci più invarianza nell’ambito delle frequenze (si veda M.G. Albanesi, S. Bertoluzza 1994).

In figura 8 si mostra un confronto per la figura Lena (uno standard per i metodi di compressione) utilizzano la base di Haar con e senza la metrica adattata al sistema visivo umano data dalla (12).

Figura 8 : raffronto tra compressione usuale e con sistema visivo

Attualmente è allo studio un metodo che fornisca l’accoppiamento tra questo approccio ed il metodo di compressione descritto sopra per la memorizzazione delle impronte digitali : risulterà essenziale avere una migliore comprensione delle caratteristiche specifiche di questa classe particolare di immagini.

Gli schemi di compressione delle immagini si basano usualmente, come abbiamo già detto, sulla presenza di informazioni rindondanti nei dati che rappreentano l’immagine stessa. Nel caso degli algoritmi di compressione frattale si ipotizza, ipotesi empirica, che ci siano regioni dell’immagine in qualche senso ‘‘simili’’ ad altre parti dell’immagine ma su scala differente.

Inizialmente le tecniche frattali sono state introdotte in computer graphics per modellizzare fenomeni naturali (montagne, alberi, foglie, conchiglie, ...). In seguito, partendo dalla teoria IFS (Iterated Function Systems) proposta da J. Hutchinson, M. Barnsley ed il suo gruppo di lavoro del Georgia Institute of Technology, hanno approfondito la possibilità della codifica di intere immagini sfruttandone l’eventuale auto-similarità (propria delle figure frattali). Dal punto di vista computazionale sono stati poi determinanti i contributi di A. Jacquin e Y. Fisher. Introduciamo le idee fondamentali. In generale sia dato uno spazio metrico completo M con distanza d(x,y), x,y Î M (il lettore disorientato pensi al piano con la distanza euclidea). Un IFS è un insieme finito di applicazioni { w_i} , w_i : S à S i=1,...,n, contrattive cioè tali che d(w_i(x),w_i(y))£ s_i d(x,y) per ogni x,yÎ S con 0 £ s_i £ 1. Consideriamo l’insieme B sottoinsiemi chiusi e limitati BÌ S dello spazio metrico considerato e la trasformazione W : B à B , W(B)= È _i w_i(B). La trasformazione W risulta essere contrattiva per B e rispetto ad una opportuna distanza definita per elementi di B . Dal Teorema delle contrazioni si deduce allora che esiste un elemento AÎ B tale che

W(A)=A=È _i w_i(A).

Tale elemento A si chiama anche insieme attrattore e può essere equivalentemente caratterizzato dal limite A=lim _{n® 8} W⁽ⁿ⁾ (B), dove l’insieme B è un qualsiasi elemento di partenza scelto in B e dove con W⁽ⁿ⁾ abbiamo indicato la composizione n-esima della trasformazione W con se stessa : W⁽⁰⁾ =W, W⁽ⁿ⁾ =W(W^(n-1)). Vediamo un semplice esempio in cui si utilizzino tre funzioni w_i che risultanmo essere tre trasformazioni affini del piano e come spazio metrico il piano euclideo.

Indichiamo con (x₁,x₂) e (y₁,y₂) due punti del piano e (x’₁,x’₂), (y’₁,y’₂) i punti trasformati, abbiamo

w₁  data da x’₁=0.5 x₁, y’₁=0.5 y₁ ; w₂  data da x’₁=0.5 x₁+0.5, y’₁=0.5 y₁ ; w₃  data da x’₁=0.5 x₁+0.25, y’₁=0.5 y₁+0.5 . L’insieme B di partenza è un quadrato e si mostrano in figura 9 alcune iterazioni successive dell’operatore W : R ² à R ², l’insieme attratore A è noto come triangolo di Sierpinski.

Barnsley ha suggerito di memorizzare al posto dell’immagine (nel caso di figura 9 il triangolo di Sierpinski) la collezione di trasformazioni utilizzate. Per riottenere l’immagine originale basta applicare la trasformazione memorizzata in modo iterativo partendo da un qualsiasi insieme per riottenere l’attrattore (o una sua approssimazione). In pratica bastano 10-12 iterazioni per avere una buona ricostruzione dell’immagine originale. Dal punto di vista matematico occorre risolvere un problema inverso : dato un certo sottoinsieme determinare l’insieme di applicazioni il cui attrattore lo approssima bene. Tale problema risulta, in generale, di ‘‘difficile’’ soluzione e rientra nell’ambito dei problemi mal posti. Mentre un IFS produce un unico attrattore, un dato attrattore può essere generato da più IFS differenti.

Le immagini naturali non sono generalmente un attrattore per un IFS, mentre singole porzioni dell’immagine sono spesso simili a parti differenti e tale relazione può essere codificata tramite un opportuno IFS. Un algoritmo di compressione frattale deve essere in grado di trovare questa similarità locale e codificarlo opportunamente. Si partiziona quindi l’immagine in blocchi R_i , i=1,..,N e blocchi più grandi D _ij, j=1,...,M (per esempio R_i siano blocchi di 4´ 4 o 8´ 8 pixels mentre D _j siano i blocchi di 8´ 8 o 16´ 16 pixels). Per ogni R_i si ricerca il dominio D tra i candidati D _j e le trasformazioni w_k (solitamente affini) che formino la trasformazione W, in modo tale da minimizzare la distanza tra R_i ´ [a,b] e w_k(D). Qui l’intervallo [a,b] è stato aggiunto per poter utilizzare anche il contrasto e la luminosità del blocco (per esempio un R_i potrebbe essere la versione ruotata, traslata e con minore luminosità di un blocco D). Il punto fondamentale, dopo aver deciso la metrica utilizzata per valutare la distanza tra immagini, consiste proprio nella ricerca dei blocchi R, D e nel trovare le trasformazioni. Date le trasformazioni ed i blocchi si codificano queste informazioni ottenendo buoni livelli di compressione.

In figura 10 si mostra un esempio di immagine ottenuta tramite questo procedimento con blocchi R e D di dimensioni citate sopra con compressione 1 :6 e SNR=27.7 (l’immagine originale è l’impronta digitale di figura 1 compresa nella galleria di immagini Matlab).

Da un confronto con l’immagine originale, figura 1, si nota la buona qualità dell’immagine ricostruita. Per la ricostruzione si utilizza il metodo iterativo con ora, però, più sistemi D,w_k. Si osservi che questa trasformata è molto assimetrica, infatti il costo computazionale è molto più elevato in fase di compressione, problemi inversi da risolvere, che in fase di decompressione, semplici iterazioni di trasformazioni affini. Contrariamente alla trasformata wavelet risulta difficile trovare indicatori di regolarità locale, come veniva fornito dai dettagli in una trasformazione multiscala, o che caratterizzino pattern differenti. In ogni caso dagli esperimenti fatti risulta che è possibile ottenere alti fattori di compressione con una buona qualità dell’immagine compressa tramite la metodologia frattale.

In questo lavoro abbiamo trattato la possibilità di utilizzare le wavelets o la metodologia frattale per la compressione di immagini raffiguranti impronte digitali.

Da una parte abbiamo la versatilità e la rapidità della trasformata wavelet discreta, mentre invece per la compressione frattale si è sperimentata la capacità di avere alti fattori di compressione mantenendo l’errore ‘‘piccolo’’. Vediamo alcuni possibili sviluppi :

• approfondire lo studio e la ricerca dell’accoppiamento tra wavelets e sistema visivo umano per migliorare la qualità dell’immagine ricostruita ed eventualmente aumentare il rapporto di compressione tramite una opportuna selezione dei dettagli che tenga conto della funzione di trasferimento dell’occhio e della particolare classe di immagini utilizzata.
• Accoppiare la trasformata wavelet con la compressione frattale sfruttando la possibilità di dedurre dai blocchi di dettagli eventuali similitudini tra parti dell’immagine, In questo modo si potrebbero sfruttare le buone proprietà di entrambi gli approcci.
• Tramite le basi di funzioni wavelets è possibile estrarre informazioni locali relative ad una determinata funzione. Ai fini dell’identificazione e della classificazione delle impronte digitali occorre approfondire quali caratteristiche devono essere estratte e come tradurre tali caratteristiche nei coefficienti dello sviluppo in una base di wavelets. La classificazione e l’identificazione potrebbero avvenire tramite l’utilizzo di una rete neuronale a più strati che permetta la combinazione delle caratteristiche morfologiche locali.
• Occorre prevedere la possibilità di garantire una maggiore segretezza ai dati utilizzati (visto l’utilizzo in ambiente investigativo) per permettere anche una trasmissione sicura degli stessi. Si dovrà quindi approntare un apposito tool di crittografia per una codifica sicura delle informazioni.

Per quanto riguarda l’implementazione dei metodi notiamo che per la parte riguardante la trasformata wavelet si è utilizzato l’apposito Wavelet-Toolbox con una notevole semplificazione. Per la compressione frattale invece non abbiamo ancora una implementazione totalmente integrata nell’ambiente Matlab a causa di una maggiore ‘‘sofisticazione’’ nelle strutture dati. Pensiamo comunque di poter produrre nel prossimo futuro un opportuno Toolbox .

esistono anche dei bollettini elettronici periodici a cui è possibile abbonarsi. Storicamente il primo, e più completo, e’ il wavelet digest edito da Wim Sweldens dei Bell Laboratories. Esistono anche siti più specialistici, per esempio per l’analisi delle tessiture (utile anche per la classificazione futura delle impronte digitali) si può ‘navigare’ fino a :

Per quanto riguarda la parte teorica sicuramente i testi di Meyer sono i più completi, ma non di facile lettura