gestione immagini - Prof. Vito Bevilacqua

POLITECNICO DI BARI
CORSO DI LAUREA
ING. ELETTRONICA SPECIALISTICA S.E.M.A.
Corso di Informatica Medica
Docente: Prof. V.Bevilacqua
GESTIONE IMMAGINI
BENEDETTO Anna Cinzia
ANNO ACCADEMICO 2008-2009
1
INDICE
1.Introduzione
1.1.Tipologie di immagini
1.2.Sequenze di immagini
2.Elaborazione immagini
2.1.Informazioni su un file di immagine
2.2.Leggere l’immagine
2.3.Scrivere dati su file di immagine
2.4.Convertire formato immagini
2.5.Visualizzare immagini
2.6.Image registration
3.Immagini DICOM
3.1.Il formato DICOM
3.2.Immagini DICOM
3.3.Immagini DICOM in MATLAB
4.Trasformate
4.1.Trasformata di Fourier
4.1.1.Visualizzare l’immagine trasformata in MATLAB
4.2.Trasformata di Fourier discreta (DFT)
4.2.1.Trasformata 2-D
4.2.2.Proprietà trasformata 2-D
4.3.DFT applicata alle immagini
4.3.1.FFT 2D
4.3.2.Analisi di uno spettro
4.3.2.1.Visualizzare uno spettro
4.3.2.2.Interpretare uno spettro
4.3.3.Applicare la FFT alle immagini
4.4.Applicazioni in MATLAB
4.5.Altri tipi di trasformata
4.5.1.Trasformata coseno discreta
4.5.2.Trasformata di Radon
5.Filtraggio
5.1.Dominio dello spazio
5.2.Dominio delle frequenze
6.Operazioni morfologiche
p.4
p.5
p.6
p.7
p.7
p.7
p.7
p.8
p.8
p.9
p.11
p.11
p.12
p.15
p.17
p.18
p.19
p.20
p.21
p.22
p.22
p.22
p.23
p.23
p.25
p.27
p.29
p.31
p.31
p.34
p.37
p.37
p.40
p.44
2
6.1.Dilatazione ed erosione
6.2.Individuare e misurare oggetti all’interno di una immagine
6.3.Selezionare un singolo oggetto in una immagine binaria
6.4.Operazioni tramite uso di Look up Table
7.Analizzare le immagini
7.1.Acquisire informazioni sul valore dei pixel e sulla
statistica delle immagini
7.2.Eliminare il rumore da una immagine
8.Utilizzare una regione di interesse (ROI)
8.1.ROI
8.2.Filtrare una ROI
p.44
p.44
p.46
p.46
p.47
p.47
p.47
p.48
p.48
p.48
3
1.INTRODUZIONE
Prima di introdurre le funzioni fondamentali per l’utilizzo di immagini in matlab, é
necessario conoscere come matlab gestisce i dati relativi alle immagini per poter operare
con questo strumento.
Matlab gestisce qualsiasi dato in forma di array, quindi anche una immagine verrà tradotta
in forma matriciale, ed ogni elemento della matrice rappresenta un pixel dell’immagine (per
le immagini a colori la matrice avrà tre dimensioni solo per la componente colore, una per
ognuno dei tre colori principali (rosso, verde e blu), oltre che le dimensioni proprie
dell’immagine).
Per individuare un punto dell’immagine, quindi, si può utilizzare un sistema discreto di
coordinate a pixel: si avrà una localizzazione del punto attraverso l’indicazione della riga e
della colonna come in una qualsiasi matrice di dati, ad esempio I(r,c) individua il valore del
pixel alla riga r colonna c.
E’ però possibile ricorrere ad un sistema di coordinate continuo per la localizzazione di un
pixel, quello spaziale, che risulta coincidente con quello a pixel per valori interi (le
coordinate spaziali del centro di un pixel coincidono con le coordinate a pixel per quel
pixel), ma che differisce nella denominazione del primo pixel ((1,1) per quello a pixel,
(0.5,0.5) per quello spaziale) e nella notazione usata per individuare le componenti
orizzontali e verticali che risultano invertite tra un sistema e l’altro( I(x,y)=I(c,r)).
4
Di default i due sistemi di coordinate coincidono e la maggior parte delle funzioni
richiedono l’immissione di dati in coordinate a pixel, ma è possibile modificare questa
informazione specificando il range dimensionale dell’immagine in due vettori.
1.1.Tipologie di immagini
::Binary images (BW): immagini bilivello nelle quali ogni pixel può assumere solo
due valori discreti: 0 (nero) e 1 (bianco)
::Indexed images (X): immagini che consistono in una prima matrice (X) i cui valori
indicano i singoli pixel e che rimandano ad una seconda matrice dei colori (map) le cui
colonne indicano i livelli di rosso, verde e blu del valore di X corrispondente al singolo
pixel.
::Grayscale images (I): le matrici corrispondenti hanno un elemento per pixel e il
valore degli elementi varia all’interno di un range di rappresentazione di grigi.
5
::Truecolor images (RGB): immagini a colori i cui pixel sono rappresentati
attraverso indicazioni sui tre colori principali – rosso, verde e blu- all’interno di un array di
tre dimensioni.
1.2.Sequenze di immagini
A volte ci si trova ad aver a che fare con sequenze di immagini tra loro collegate, come per
esempio slice di risonanze magnetiche. Per gestire le informazioni fornite matlab utilizza
array multidimensionali: ad esempio trasforma una sequenza di p immagini di dimensione a
x b in un array axbxp.
6
Molte funzioni operano su array multidimensionali, ma alcune non sono in grado di
riconoscere in quel tipo di dato una sequenza di immagini, quindi è necessario usare una
particolare sintassi e rispettare alcuni limiti di utilizzo (vedere “Working with Image
Sequences”).
2.ELABORAZIONE IMMAGINI
E’ buona norma come prima operazione da compiere, nel caso si voglia elaborare una
immagine in matlab, ripulire il workspace dalle variabili precedentemente utilizzate e
chiudere le eventuali immagini aperte tramite il comando close all.
2.1.Informazioni su un file di immagine
Per conoscere le informazioni riguardanti un file di immagine come il tipo, la data di
creazione, le dimensioni in termini di pixel e byte, il numero di immagini (se si sta
analizzando una sequenza di immagini), si utilizza
imfinfo(‘nomefile.estensione’)
2.2.Leggere l’immagine
Per leggere una immagine, di qualsiasi estensione, e tradurla in una matrice di dati I si usa
I=imread(‘nomefile.estensione’).
Nel caso in cui il file si composto da più immagini, questa funzione carica informazioni solo
dalla prima. Ad esempio per leggere le informazioni di 27 immagini di un file .tiff e
caricarle in una matrice di 4 dimensioni:
mri = zeros([128 128 1 27],'uint8');
for frame=1:27
[mri(:,:,:,frame),map] = imread('mri.tif',frame);
end
2.3.Scrivere dati su file di immagine
Dopo aver tradotto una immagine come matrice di dati, e magari averla modificata, per
poter osservare l’immagine relativa bisognerà compiere una operazione inversa rispetto la
precedente attraverso la funzione
imwrite(I,’nomefileimmagine.estensione’)
7
In base alla estensione scelta si ha la possibilità di aggiungere ulteriori parametri come la
qualità e compressione per il formato jpeg (vedere pagina di riferimento della funzione
imread per i particolari).
2.4.Convertire formato immagini
Per convertire il formato di una immagine è necessario caricarla come file di dati (attraverso
imread) e successivamente utilizzare imwrite per salvarla con una estensione diversa. Un
esempio è dato dalle seguenti linee di codice:
moon_tiff = imread('moon.tif');
imwrite(moon_tiff,'moon.jpg')
2.5.Visualizzare immagini
La visualizzazione dei file di immagini è affidata a due funzioni imshow e imtool.
La funzione imshow è utilizzata per visualizzare immagini di tutti i tipi supportati da matlab
e può essere espressa in due forme, la prima utilizza la funzione imread
immagine = imread('nomefile.tif');
imshow(immagine);
e l’altra richiama direttamente il nome del file
imshow(‘nomefile.tif')
Utilizzando la seconda espressione, però, MATLAB si limita a far visualizzare l’immagine
ma non ne importa i dati nel workspace. A questo scopo è necessario utilizzare le funzione
che assegna i dati dell’immagine alla variabile “immagine” se si vuole avere la possibilità di
richiamare l’immagine in seguito
immagine= getimage
La funzione imshow apre l’immagine di default al 100% di zoom, la racchiude in margini e,
in caso di utilizzo consecutivo con immagini diverse, le immagini si rimpiazzano, ma tutte
queste impostazioni sono modificabili.
La funzione imtool è utilizzata per visualizzare le immagini, ma apre anche una finestra
attraverso la quale poter analizzare molte delle caratteristiche dell’immagine, quali la
posizione di un pixel, il dato matriciale corrispondente a un gruppo di pixel, la distanza tra
pixel, visualizzare informazioni base sull’immagine (estensione e dimensioni, oltre a
informazioni su ciò che è rappresentato) e modificarla. Anche questa funzione, come la
precedente può essere usata in due diverse forme
immagine = imread('nomefile.tif');
imtool(immagine);
e
imtool('nomefile.tif')
8
Anche in questo caso la seconda forma non permette l’importazione dei dati dell’immagine
nel workspace di MATLAB.
Questa funzione, aprendo una finestra di dialogo che permette di modificare l’immagine,
permette di modificare la visualizzazione, importare ed esportare dati da e al workspace, di
salvare i dati, e compiere molte altre attività.
L’espressione di utilizzo di queste funzioni varia anche se di poco da tipo a tipo di
immagine, perché si deve considerare la diversa struttura dei dati associati alle immagini nei
vari formati.
2.6.Image registration
Si tratta di un processo attraverso il quale è possibile allineare immagini diverse raffiguranti
uno stesso oggetto in modo da ricostruirlo più chiaramente, paragonando le caratteristiche
delle singole immagini.
Il processo consiste nell’allineare all’immagine base le immagini di input attraverso
modifiche spaziali (come rotazioni, estrazione di porzioni di immagine e così via) in modo
da renderle confrontabili.
Uno strumento usato per questo tipo di confronto è l’uso di punti di mappatura: per allineare
due immagini si scelgono alcuni punti che individuano caratteristiche uguali o anche i
margini delle due immagini, dai quali si parte nell’analizzare il resto delle immagini.
Ovviamente può essere necessario iterare questa operazione per ottenere buoni risultati.
I passi attraverso i quali si sviluppa questo processo sono i seguenti:
9
Nel caso si voglia fare lo stesso tipo di registrazione per più di due immagini, si può
automatizzare il processo creando uno script in modo che, fornendo l’immagine base e
quelle di input e lanciando il tool di selezione dei punti di mappatura, l’operazione diventi
immediata. Le operazioni necessarie all’operazione sono compiute attraverso l’uso della
funzione cpselect.
moon_base = imread('moon.tif');
moon_input = moon_base;
cpselect(moon_input, moon_base);
(L’uso di questa funzione si accompagna a quello dell’opzione wait, che blocca le
operazioni di MATLAB finché è in corso il controllo dei punti di mappatura, evitando che il
controllo sui punti di mappatura scelti sia parziale). Lo script creerà una TFORM, una
struttura di MATLAB che contiene tutti i parametri necessari ad ottenere delle
trasformazioni spaziali 2-d (traslazioni, scalamento delle dimensioni, rotazioni..) e la
fornirà, insieme con l’immagine di input, alla funzione imtransform .
Si potrebbe migliorare la scelta dei punti di mappatura (che fino ad ora si è ipotizzato essere
fatta sulla base della semplice vista) attraverso l’uso della cross-correlazione.
Per utilizzare la cross-correlazione, si passano i punti di mappatura dell’immagine base e di
quella di input e le immagini stesse alla funzione cpcorr.
10
input_pts_adj= cpcorr( input_points, base_points, input, base).
Questa funzione definisce un array 11 per 11 attorno ad ogni punto di controllo
dell’immagine di input e di ogni punto di mappatura dell’immagine base e calcola la
correlazione tra i pixel delle regioni così individuate. In seguito, in base al valore di
correlazione, viene ottimizzata la scelta dei punti di controllo attraverso una riduzione del
loro numero.
3.IMMAGINI DICOM
3.1.Il formato DICOM
Ultimamente si sta presentando sempre più spesso la necessità di comunicare informazioni
di tipo medico. La possibilità di scambiare immagini mediche corredate da varie
informazioni attraverso reti telematiche, però, fa sorgere l'esigenza di uno standard di
comunicazione globale, in modo tale da offrire al radiologo la possibilità di consultare tutti
gli esami dello stesso paziente contemporaneamente ed in un solo posto, ed al medico di
reparto di accedere immediatamente alle immagini, attraverso un collegamento in rete
locale. Uno standard, infatti, consiste in un insieme di norme fissate allo scopo di ottenere
l'unificazione delle caratteristiche di una determinata prestazione o processo tecnologico, da
chiunque o comunque prodotto. Il formato DICOM (Digital Imaging and Communications
in Medicine) rappresenta un modello di tale standardizzazione.
Fin dalla sua prima apparizione questo standard ha rappresentato un’importante novità. Da
un lato, ha definitivamente stabilito la non divisibilità tra i dati dei pixel dell’immagine
medica e i dati descriventi il procedimento che ha portato alla formazione dell’immagine
stessa. A tal fine, lo standard DICOM include entrambe le informazioni in un unico file che
non rappresenta dunque un’immagine così come potrebbe essere descritta da altri formati,
ma al tempo stesso un’immagine, a chi e a cosa quest’immagine fa riferimento ed in che
modo è stata ottenuta. Dall’altro, lo standard ha promosso l’interconnessione tra le
apparecchiature di diagnostica medica: non più apparecchiature di diagnostica a sé stanti,
ma apparecchiature di diagnostica, e in alcuni casi di terapia, postazioni di visualizzazione e
di refertazione che possono essere tra loro collegate e che sono in grado di condividere
immagini e dati correlati, dispositivi di archiviazione e stampanti.
L’interoperabilità tra le apparecchiature è garantita in quanto lo standard DICOM definisce
un protocollo di comunicazione, ovvero un insieme di regole da rispettare nel momento in
cui le apparecchiature dialogano tra loro, che è basato sul protocollo di comunicazione della
rete internet. L’adozione del protocollo DICOM è il solo requisito richiesto affinché
un’apparecchiatura possa far parte della rete. Il modello di comunicazione è quello
Client/Server che lo standard ridefinisce User/Provider. Sono “provider” i computer che
sono in grado di fornire un servizio, mentre sono “user” quelli che richiedono un
11
determinato servizio. I principali “servizi”, ovvero le possibili azioni che è possibile
intraprendere su una rete DICOM, sono di seguito elencati:
- Storage: è il servizio richiesto da un’apparecchiatura per “archiviare” immagini; è
utilizzato tipicamente da una modalità di acquisizione per inviare immagini ad un server di
archiviazione .
- Storage Commitment: è un servizio di storage con in più un messaggio di “avvenuta
archiviazione” inviato dal Provider del servizio di storage allo User, cosicché questo possa
“cancellare” localmente le immagini.
-Query/Retrieve: è il servizio richiesto quando da un’apparecchiatura se ne interroga
un’altra al fine di conoscere la lista delle immagini su questa disponibili ed eventualmente
richiamarle.
- Modality Worklist: è il servizio in grado di gestire la lista degli esami da acquisire per
ciascun paziente. Ogni esame della lista viene programmato e il suo completamento viene
reso noto al sistema che provvede ad aggiornare i dati. E’ possibile solo in dipartimenti
dotati di un sistema di prenotazione/accettazione informatizzato integrato nella rete
DICOM.
- Print: per la stampa delle immagini da una modalità di acquisizione o da una stazione di
visualizzazione. Non è detto che un’apparecchiatura DICOM-compliant disponga di tutti i
servizi sopra elencati. La disponibilità’ o meno di alcune funzioni deve essere verificata
all’atto dell’acquisto. Il costruttore dell’apparecchiatura è tenuto a riportare in un
documento di conformità (Conformance Statement) i dettagli dell’implementazione.
3.2.Immagini DICOM
I dati radiologici rappresentabili come immagini o le immagini vere e proprie che vengono
archiviate secondo lo standard DICOM sotto forma di file vengono comunemente chiamate
immagini DICOM. L'errore più comune che viene fatto nell'interpretazione del termine è
che queste siano assimilabili ad altri formati di compressione dell’immagine come JPEG o
GIF. In verità lo standard DICOM applicato alla codifica dei file non è nient’altro che un
metodo per incapsulare i dati e per definire come questi debbano essere codificati o
interpretati, ma non definisce alcun nuovo algoritmo di compressione. La maggior parte
delle volte, l'immagine viene archiviata in forma non compressa, secondo la codifica con la
quale viene prodotta, ma esistono molti software che sono in grado di produrre o
interpretare file DICOM contenenti dati compressi secondo vari algoritmi.
Le informazioni che il DICOM registra all’interno del file che rappresenta un’immagine, in
aggiunta a quelle relative ai pixel, sono sistemate in gruppi all’interno dell’header, secondo
il seguente ordine: Paziente, Studio, Serie e Immagine.
Nel gruppo Paziente sono riportate le informazioni anagrafiche della persona sottoposta a
indagine medica.
12
Nel gruppo Studio sono definite le caratteristiche dell’indagine diagnostica che il paziente
deve effettuare e che può richiedere esami di differenti modalità.
Nel gruppo Serie è descritta ciascuna collezione di immagini proveniente da una singola
modalità diagnostica ed i relativi parametri di acquisizione.
13
Infine, nel gruppo Immagine sono definiti gli attributi dei pixel che compongono
l’immagine, quali dimensione della matrice, profondità del pixel, rappresentazione del pixel
e interpretazione fotometrica.
L’immagine seguente rappresenta un esempio di una tipica immagine DICOM :i primi 794
bytes sono utilizzati per l’header ed i successivi 19838 byte per l’immagine in sé. La
dimensione dell’header non è fissa, ma dipende da quante informazioni si vogliono
immagazzinare.
L’immagine seguente mostra una descrizione più dettagliata dell’header: è immediatamente
indicata una caratteristica di questo tipo di file, una sorta di identificativo, e cioè che
DICOM richiede un preambolo di 128-byte, di solito sono tutti impostati allo zero, seguita
dalle lettere “D”, “I”, “C”, “M”.
Questo preambolo è seguito dalle informazioni dell’header, che sono organizzate in gruppi.
14
La quantità di informazioni da inserire in una immagine DICOM, dipende dal tipo di
immagine che si sta trattando: ad esempio l’immagine che si sta analizzando è del tipo
“MR”, quindi deve contenere gli elementi necessari a descrivere il processo temporale di
risonanza; la mancanza di tali informazioni in questa immagine è una violazione dello
standard DICOM.
3.3.Immagini DICOM in MATLAB
Le immagini DICOM sono immagini di tipo medico, nelle quali sono contenute
informazioni di varia natura sul paziente e sulla rilevazione dell’immagine, oltre
all’immagine in sé. Per poter leggere questo tipo di informazioni si procede con la seguente
riga di codice
info=dicominfo(‘nomefile.dcm’).
Si otterrà un elenco di informazioni come in seguito:
info =
Filename: [1x47 char]
FileModDate: '24-Dec-2000 19:54:47'
FileSize: 525436
Format: 'DICOM'
FormatVersion: 3
Width: 512
Height: 512
BitDepth: 16
ColorType: 'grayscale'
15
SelectedFrames: []
FileStruct: [1x1 struct]
StartOfPixelData: 1140
MetaElementGroupLength: 192
FileMetaInformationVersion: [2x1 double]
MediaStorageSOPClassUID: '1.2.840.10008.5.1.4.1.1.7'
MediaStorageSOPInstanceUID: [1x50 char]
TransferSyntaxUID: '1.2.840.10008.1.2'
ImplementationClassUID: '1.2.840.113619.6.5'
.
.
.
Alle volte, oltre a questo tipo di informazioni, ci sono anche altre informazioni di tipo
privato, che risultano quindi inaccessibili e sono restituite dalla finzione dicominfo come un
generico nome sulla base delle informazioni sul gruppo di appartenenza di quel dato: ad
esempio se il dato è contenuto nel gruppo 0009 all’elemento 0006, dicominfo creerà un
nome del tipo “Private_0009_0006”.
Tutti campi di informazione che è possibile specificare in immagini DICOM sono catalogati
all’interno di una libreria di Matlab; nel caso in cui si voglia inserire dei campi che non
compaiono in questa libreria, è possibile aggiornare la libreria creando una propria versione
modificabile.
Per leggere una immagine DICOM in forma di matrice di dati si possono usare due modi: il
primo fornendo il nome del file che si vuole leggere
I=dicomread(‘nomefile.dcm’)
il secondo attraverso le informazioni aggiuntive sul paziente e sulla rilevazione
dell’immagine
info = dicominfo(‘nomefile.dcm’);
I = dicomread(info).
Per visualizzare immagini DICOM in scala di grigi si utilizzano le funzioni imshow e
imtool, ma è necessario in entrambi i casi fare preventivamente un’autoscale dell’immagine
che, essendo l’immagine DICOM a 16 bit, non potrebbe essere visualizzata altrimenti.
imshow(I,'DisplayRange',[ ])
Per scrivere i dati dell’immagine in un file di immagine DICOM utilizzo la funzione
dicomwrite(I,'nomeimmaginedicom.dcm')
che scrive i dati dell’immagine I nel file DICOM specificato.
Una volta chiamata questa funzione, essa stessa di default include nell’immagine DICOM in
creazione il minimo numero di informazioni richieste dal tipo di oggetto DICOM (IOD) che
si sta creando; è possibile creare come oggetti una “Secondary Capture” (di default),
“risonanza magnetica ” e “tomografia”. Altra possibilità che si ha per definire i campi di una
16
immagine DICOM è quella di utilizzare la struttura di dati di un file già esistente e passarla
alla funzione dicomwrite. Di seguito sono descritte le righe di codice necessarie per questo
tipo di operazione: di passa alla funzione in esame l’insieme di informazioni che sono state
ottenute attraverso dicominfo.
info = dicominfo('CT-MONO2-16-ankle.dcm');
I = dicomread(info);
dicomwrite(I,'ankle.dcm',info)
Il formato DICOM, come si è già detto, prevede una struttura nella quale vengono allocate
prima le informazioni dell’header e poi i dati dell’immagine. E’ possibile includere un
“valore di rappresentazione” VR nell’header, oppure lasciare che sia settato
automaticamente sulla base dei dati della libreria: per specificare il VR è necessario inserire
il parametro ‘explicit’ nel momento in cui viene chiamata la funzione dicomwrite.
L’immagine seguente mostra l’header con e senza il VR specificato.
Dato che le immagini di tipo DICOM contengono anche informazioni riservate sul paziente,
può essere utile rendere anonimi i dati a disposizione andando a modificare le informazioni
proprie del file tramite la funzione
dicomanon(immagine_in,immagine_out)
che elimina le informazioni riservate sul paziente dal file immagine_in e crea un nuovo file
(immagine_out) privo di queste informazioni, e quindi divulgabile.
4.TRASFORMATE
In genere la rappresentazione matematica di un immagine è una funzione di due variabili
spaziali f(x,y). Il valore della funzione in un determinato punto rappresenta l'intensità
dell’immagine in quel punto: questo rappresenta il dominio spaziale.
Parlando di trasformata di Fourier, invece, ci si riferisce ad una rappresentazione
matematica alternativa e che fornisce uno strumento molto potente per analizzare ed
elaborare le immagini, ma per poter “capire” il contenuto occorre guardare la trasformata di
una immagine in modo del tutto diverso da come si “guarda” una immagine. La trasformata
di Fourier è una rappresentazione di un'immagine come una somma di esponenziali
17
complessi di diversa ampiezza, frequenza e fase, e questo rappresenta il dominio delle
frequenze, le frequenze spaziali, in questo modo ogni pixel contiene contemporaneamente
un'informazione relativa all'ampiezza (cioè al fatto che l’immagine contenga o meno una
struttura periodica) e una relativa alla fase (cioè dove le eventuali strutture periodiche sono
localizzate). La frequenza spaziale fornisce una rappresentazione legata al livello di
dettaglio presente in una immagine; matematicamente, può essere vista come un profilo di
luminanza sinusoidale variabile fra un massimo (bianco) ed un minimo (nero): un'immagine
sfocata (quindi con basso dettaglio) può essere considerata a bassa frequenza, mentre
un'immagine che contenga solo profili può essere considerata ad alta frequenza spaziale.
z=f(x,y)=A sin(2π (ux+vy)+φ)
Nei processi di elaborazione delle immagini si può decidere di utilizzare sia il dominio
spaziale che quello frequenziale, in quanto esistono algoritmi specifici per questo compito
in entrambi i domini; spesso però si opta per un approccio frequenziale che risulta meno
oneroso e molto più veloce, come accade per la FFT (Fast Fourier Transform).
La trasformata, quindi risulta utile per una vasta gamma di applicazioni di elaborazione
delle immagini, come il miglioramento della qualità, l’analisi e la compressione.
4.1.Trasformata di Fourier
Presa la funzione f(m,n) di due variabili spaziali discrete m ed n, la sua trasformata di
Fourier bidimensionale, che definisce la rappresentazione nel dominio delle frequenze di
f(m,n), è definita da
F (w1, w2) = ∑∑ f (m, n)e− jw1me− jw2n
m
n
18
dove w1 e w2 sono le variabili in frequenza e la F(w1,w2) risulta una funzione complessa e
2π-periodica. Si può notare che F(0,0) rappresenta la somma delle componenti di f(m,n) ed
è chiamata componente costante della trasformata ma, dato che gli indici di MATLAB
partono da 1 e non da 0, il programma riconoscerà F(1,1) come componente costante.
Una volta ottenuta la trasformata di una immagine, e magari elaborata con filtri,è possibile
ricostruire l’immagine nel dominio spaziale attraverso l’antitrasfomata di Fourier, definita
come
π
π
1
f (m, n) = 2 ∫ ∫ F (w1, w2)e jw1me jw2ndw1dw2
4π w1 =−π w2 =−π
Questa forma indica che è possibile rappresentare f(m,n) come somma di infiniti
esponenziali complessi di diversa frequenza.
4.1.1.VISUALIZZARE L’IMMAGINE TRASFORMATA IN MATLAB
Considero la funzione spaziale f(m,n) nella figura seguente che descrive, anche se in modo
continuo e non discreto per facilità, un rettangolo di valore 1 all’interno e 0 all’esterno.
Applicando la trasformata di Fourier a questa immagine, sarà possibile visualizzare il
modulo della funzione trasformata |F(w1,w2)|.
Si può notare come il processo di trasformata
19
E’ facilmente individuabile la F(0,0) nel picco centrale (somma dei valori della f(m,n)); si
osserva anche come le alte frequenze orizzontali presentino una ampiezza sicuramente
maggiore rispetto le alte frequenze verticali, come risultato del fatto che la f(m,n) si
sviluppa maggiormente lungo la direzione verticale che lungo quella orizzontale. Altro
modo per visualizzare il modulo della la componente frequenziale di una immagine è la
visualizzazione logaritmica log|F(w1,w2)|.
Questo tipo di visualizzazione aiuta a mettere in evidenza i dettagli circa le ampiezze
dell’immagine trasformata lì dove sono molto vicine allo zero.
E’ preferibile utilizzare questo tipo di visualizzazione, altrimenti la componente continua, di
ampiezza molto maggiore al resto delle componenti,oscurerebbe il resto dello spettro.
(Notare che per visualizzare il log devo precedentemente scrivere una riga di codice che
faccia il modulo.)
4.2.Trasformata di Fourier discreta (DFT)
Dovendo operare con la trasformata attraverso calcolatori, ci si trova ad operare con una sua
forma discreta, la DFT. Questa forma ben si presta all’utilizzo di calcolatori dato che i suoi
ingressi e uscite sono valori discreti, e inoltre questa forma di trasformazione è supportata
da un algoritmo veloce di risoluzione della trasformata, la FFT (Fast Fourier Transform).
La DFT è utilizzata per funzioni discrete definite diverse da zero solo in una regione finita
M x N dove 0≤m≤M-1 e 0≤n≤N-1.
Le relazioni della DFT, diretta su un’area M x N e inversa su una N x M, sono le seguenti
20
M −1 N −1
F( p, q) = ∑∑ f (m, n)e− j(2π / M ) pme− j(2π / N)qn
m=0 n=0
1 M −1 N −1
F ( p, q) =
F ( p, q)e j (2π / M ) pme j (2π / N )qn
∑∑
MN p=0 q=0
Si può notare come le F(p,q) risultino essere i coefficienti della f(m,n) e, inoltre, siano parte
dei coefficienti della trasformata di Fourier F(w1,w2) introdotta in precedenza.
La FFT è un algoritmo sviluppato per il calcolo veloce della DFT, grazie al quale si ha una
riduzione rilevante del numero di operazioni da compiere nel calcolo della DFT. Per il
calcolo di ciascun termine della DFT bisogna eseguire N moltiplicazioni complesse e N
somme (dato che le moltiplicazioni richiedono un tempo di calcolo molto maggiore rispetto
alle somme è possibile trascurare il tempo impiegato per il calcolo di queste ultime). Il
calcolo degli N valori della DFT richiederà un numero di operazioni pari a N2
(moltiplicazioni complesse). Tuttavia, se si opera una suddivisione dell’insieme dei
campioni in due sottoinsiemi contenenti ciascuno N/2 campioni, la DFT applicata a ciascun
sottoinsieme, in questo caso per il calcolo di tutti gli N/2, richiederà un numero di
operazioni pari a (N/2)2. Per il calcolo di tutti i campioni saranno necessarie 2*(N/2)2
operazioni. Questo processo di suddivisone può essere iterato fino ad arrivare a lavorare su
due campioni per volta, fatta l’ipotesi che il numero di campioni è potenza di due (N=2k).
Essendo N=2k al più si avrà un numero di successivi dimezzamenti pari a log2N che sono
proprio pari al numero di passi necessari per eseguire la FFT. Dato che ad ogni passo al
massimo saranno necessarie N/2 moltiplicazioni allora FFT richiederà al massimo
(N/2)*log2N contro le N2 richieste nel calcolo della DFT.
Le funzioni MALAB che implementano l’algoritmo di FFT sono fft (per DFT
monodimensionali), fft2 (per DFT bidimensionali come immagini) e fftn (per DFT
multimensionali), mentre le ifft, ifft2 e ifftn implementano la DFT inversa.
4.2.1.TRASFORMATA 2-D
Invece di aver a che fare con dati di una dimensione che indicano una variazione di
ampiezza nel tempo, si può avere a che fare con immagini che indicano variazioni di
ampiezza nello spazio, e quindi si necessita dell’equivalente 2D della DFT, che richiede il
calcolo di un numero pari alla dimensione dell’immagine di coefficienti, e quindi l’utilizzo
della FFT. Per fare ciò si divide la DFT 2D in due DFT 1D, ognuna implementata
attraverso algoritmo FFT.
21
Considerando che i coefficienti da calcolare nella trasformata sono numeri complessi, si
calcolerà una matrice delle parti reali, e una delle parti immaginarie.
Il procedimento consiste nell’eseguire la trasformata dei punti della prima riga e scriverne la
parte reale e la parte immaginarie nelle due matrici. Una volta fatto ciò per tutte le righe, si
avranno a disposizione delle “matrici intermedie” della trasformata. A questo punto si
prendono i punti della prima colonna della matrice reale e di quella immaginaria e si calcola
la FFT e il risultato è scritto nella prima colonna delle rispettive matrici. Una volta ripetuto
il tutto per le restanti colonne si è ottenuta la DFT 2D dell’immagine di partenza.
Dato che un segnale reale ha spettro di ampiezza pari e di fase dispari e che se è finito nel
tempo ha spettro frequenziale periodico e infinito, per rappresentare in frequenza una
immagine NxN basterebbero N righe e N/2 +1 colonne; quindi occorre minor spazio per
memorizzare una immagine nel dominio delle frequenze.
4.2.2.PROPRIETA’ DELLA TRASFORMATA 2-D
_Separabilità_ f(x,y)= fx(y)+fy(x)
_Linearità_ Fa corrispondere a combinazioni lineari nel dominio (x,y) combinazioni lineari
nel dominio (w1,w2)
_Simmetria complessa e coniugata_
_Scalamento in ampiezza nel dominio spaziale implica scalamento in ampiezza inverso
nel dominio frequenziale_
_Shift di posizione nel dominio spaziale implica shift di fase nel dominio frequenziale_
_Teorema di convoluzione_ La trasformata 2D del prodotto di convoluzione di due
funzioni nello spazio è uguale al prodotto delle trasformate nel dominio delle frequenze e
viceversa.
_Teorema di Parseval_ equivalenza energia nel dominio delle frequenze e dello spazio.
4.3.DFT applicata alle immagini
4.3.1.FFT 2D
La trasformata, per quanto detto fino ad ora, non è quindi uno strumento che aiuta ad
analizzare quello che una immagine mostra, anzi, non c’è immediata corrispondenza tra una
immagine e il suo spettro; quello che c’è di rilevante nell’utilizzare questo strumento è la
facilità di elaborazione e minor uso della memoria (dipende dalla struttura stessa della FFT
che procede iterativamente secondo uno schema a farfalla e ad ogni iterazione, dato che i
campioni utilizzati non servono più, si sovrappongono i risultati ai campioni stessi).
Molte sono le considerazioni di tipo matematico che sono alla base dell’algoritmo della FFT
e della complessità numerica associata al suo utilizzo che bisogna tener presente nello
22
studio dell’applicazione della FFT. Prima cosa da fare è capire quando è possibile applicare
la FFT
Per definizione di algoritmo FFT, si deve agire su immagini che forniscano in ingresso un
numero di campioni che sia una potenza del 2 (di solito si ha a che fare con immagini
256x256 oppure 4096x4096); se l’immagine non presenta queste caratteristiche, per poter
utilizzare la FFT, si deve forzare la dimensione dei dati in ingresso ad essere quella
desiderata attraverso una operazione di padding che aggiunge nuovi pixel con un livello di
grigio costante (di solito 0). E’ inoltre preferibile avere in ingresso alla trasformata 2D
matrici quadrate.
4.3.2.ANALISI DI UNO SPETTRO
4..3.2.1.Visualizzare uno spettro_
Primo problema che ci si pone è come visualizzare una trasformata che, per definizione,
altro non è che una somma di numeri complessi. Il modo più semplice è attraverso il
modulo, ma, dato che un numero complesso è ben altro che il semplice modulo, ci si
limiterà a considerarlo solo per la visualizzazione e non per i calcoli. L’informazione
associata alla fase è, infatti, essenziale per la struttura dell’immagine, indicando dove le
strutture periodiche individuate con la DFT sono collocate.
Alto problema è la scelta della scala di grigi che, come si detto prima dovrà essere una scala
non lineare, ma logaritmica, in modo tale che la componente continua dell’immagine non
oscuri il resto dell’immagine. Spesso è necessario ridurre la dinamica dello spettro in modo
da rendere possibile la visualizzazione di molti più dettagli. Per fare ciò si utilizza una
funzione del tipo D(u,v)= c log(1+F(u,v)) che permette, tramite la scelta della c di far
ricadere i valori trasformati nel range voluto. (ad esempio se si vuole la rappresentazione in
scala di grigi, si avrà che la massima ampiezza R assunta dall’immagine coincida con 255,
cui c=(L-1)/logR.
Ultimo problema è la visualizzazione MATLAB dell’immagine trasformata, infatti se si
calcola la trasformata di Fourier discreta di un’immagine le frequenze basse sono
rappresentate vicino ai quattro vertici e le frequenze alte al centro.
23
Se i quattro quadranti sono scambiati in diagonale si ottiene una visualizzazione più
intuitiva.
Questo comportamento è spiegabile considerando che la trasformata discreta DFT e la sua
inversa sono periodiche e la periodicità della DFT dà luogo ad una interessante
interpretazione geometrica. Nel caso1-D, il campione F(N) = F(0) è ovviamente contiguo a
F(N-1). I campioni possono quindi essere pensati calcolati per valori disposti non su una
linea retta ma su un cerchio, il cosiddetto “anello di Fourier”. Nelcaso2-D, la matrice
rappresentativa dell’immagine è proiettata sul cosiddetto “toro di Fourier”. Anche se la
F(w1,w2) si ripete infinite volte, solo gli MxN valori compresi in un solo periodo sono
necessari per ottenere f(x,y) da F(w1,w2). Quindi solo un periodo della trasformata è
necessario per specificare completamente F(w1,w2)nel dominio delle frequenze.
Se f(x) è reale, F(w)è inoltre dotata di simmetria coniugata, quindi F(w1,w2)=F*(-w1,-w2) e
|F(w1,w2)|=|F(-w1,-w2)|. Inoltre la proprietà di periodicità indica che F(w1) ha un periodo
paria N, e la proprietà di simmetria mostra che il modulo della DFT è centrato nell’origine,
infatti F(w)=F(w+N) e |F(w)|=|F(-w)|. Pertanto nell’intervallo [0, N-1] sono in realtà
presenti due semiperiodi della trasformata, costituiti, rispettivamente, dai campioni da 0 a
N/2, e dalla replica degli N/2 campioni presenti nel semiperiodo a sinistra dell’origine.
24
4.3.2.2.Interpretare uno spettro
Innanzitutto, bisogna dire che trasformando con Fourier si perde la localizzazione del
singolo pixel: ogni componente di Fourier ha una “sua” immagine che corrisponde
all’antitrasformata.
Nel caso di immagini, le sinusoidi di cui è composta la trasformata oltre ad essere
individuate da ampiezza, frequenza e fase, sono anche caratterizzate dalla direzione,
individuata nella posizione dei punti sulle matrici delle trasformate. (Bisogna congiungere
questi punti con quelli che individuano la continua, per individuare la direzione della
sinusoide).
Ad esempio, nell’immagine seguente, collegando i punti che individuano le sinusoidi con il
punto che individua la continua, è facile individuare la direzione della sinusoide.
Ulteriore analisi che è possibile fare è una giustificazione della posizione dei punti che
individuano le sinusoidi. All’immagine superiore di destra, ad esempio, corrisponde uno
spettro che individua due frequenze principali, oltre la continua, e i punti cadono su una
linea verticale che passa dal centro, perché la variabilità dell’immagine è massima se la si
attraversa in verticale.. Inoltre, supponendo che la massima frequenza rappresentabile è una
linea di un pixel (Fmax=1/pixel), dallo spettro in analisi si dedurrà che la frequenza dei
25
punti in esame sarà la meta della frequenza massima, e quindi le fasce orizzontali sono di
larghezza 2 pixel (F=1/2pixel=Fmax/2).
E’ importante analizzare come la trasformata operi con componenti armoniche e impulsi.
Ecco alcuni esempi di trasformata di Fourier su semplici funzioni periodiche:
(1)
(2)
(3)
Dalla figura 1, si osserva immediatamente che ad una componente armonica in (x,y)
corrisponde una coppia di impulsi simmetrici e coniugati in (w1,w2). Questo
comportamento era facilmente prevedibile, conoscendo il comportamento della trasformata
di Fourier: sulla base di questo si può affermare che un’altra caratteristica del passaggio da
un dominio all’altro può essere la corrispondenza ad una componente periodica in (w1,w2)
di due impulsi traslati in (x,y). Ulteriore caratteristica risultate dalla figura 2 è la
corrispondenza ad operazioni di modulazione in ampiezza di operazioni di traslazione nel
dominio delle frequenza e sulla base di ragionamenti generali si deduce che vale il
viceversa.
26
Queste proprietà permettono di aggirare molti degli ostacoli riguardati l’elaborazione di
immagini potendo andare ad agire sulle proprietà dell’immagine con facilità: una
applicazione dell’ultima proprietà esposta è quella della rappresentazione dello spettro di
una immagine, infatti, come già accennato, si utilizza una rappresentazione logaritmica per
non oscurare lo spettro ad opera della componente continua, e in più si preferisce porre la
continua al centro della visualizzazione per comodità attraverso una traslazione della DFT
(figura 3).
4.3.3.APPLICARE LA FFT ALLE IMMAGINI
Le proprietà evidenziate nel paragrafo precedente sono sfruttabili per il miglioramento o la
compressione delle immagini, ad esempio l’immagine seguente mostra una immagine sulla
quale sono evidenti gli effetti della luce sovrapposti, che sono stati assimilati ad un segnale
periodico, e quindi individuabili come impulso in frequenza; filtrato questo impulso
l’immagine risulta essere migliorata. Il particolare mostra appunto uno zoom dello spettro
delle due immagini.
Ci si può anche trovare nelle condizioni di dover rendere più nitida un’immagine sfocata;
ricordando che una immagine sfocata è considerabile a bassa frequenza e i contorni ad alta
frequenza, bisognerà filtrare l’immagine in modo opportuno come mostra l’immagine
seguente.
27
Oltre a queste applicazioni finalizzate al miglioramento dell’immagine in sé, la FFT è molto
sfruttata anche per eseguire convoluzioni di immagini con operatori grandi risultando la
FFT molto più veloce: per eseguire convoluzione tra una immagine 512x512 e un operatore
50x50 la FFT risulta più veloce di 20 volte.
L’implementazione MATLAB di questa applicazione della FFT è descritta nel paragrafo
dedicato.
Una ulteriore applicazione che vede l’impiego della FFT è il riconoscimento di una forma
assegnata (pattern-recognition): data una figura di partenza si vuole individuare una zona al
suo interno che si sa come è fatta.
La tecnica che si utilizza è quella della correlazione, implementata come convoluzione tra
l’immagine di partenza e l’immagine target da individuare. Si usa la correlazione perché
questa fornisce una stima della somiglianza tra due segnali.
Ipotizzando di voler individuare all’interno di un testo in una immagine 256x256 la lettera
“a” di 14x11, prima di tutto bisogna trasformare l’immagine target in un operatore
funzionale ruotandola di 180 gradi.
28
Ora bisogna rendere possibile la convoluzione rendendo le immagini della stessa grandezza
tramite lo zero-padding e, visto che si opererà con FFT, questa grandezza dovrà essere una
potenza del 2.
L’immagine convoluta Y, data dal prodotto delle matrici reali e immaginarie prodotte dalla
DFT, è composta da pixel la cui intensità è proporzionale a quanto l’immagine target H sia
simile a quella in ingresso X, quindi per individuare la figura cercata basterà individuare il
pixel più luminoso (toni del rosso).
Anche l’implementazione MATLAB di questo tipo di applicazione della FFT è descritta nel
paragrafo seguente.
4.4.Applicazioni in MATLAB
Un primo esempio di codice MATLAB che si può fornire è quello relativo alla semplice
applicazione di FFT e visualizzazione dello spettro.
1) Costruzione di una semplice immagine
f = zeros(30,30);
f(5:24,13:17) = 1;
imshow(f,'InitialMagnification','fit')
2) Visualizzazione imagine trasformata tramite la funzione fft2.
F = fft2(f);
F2 = log(abs(F));
imshow(F2,[-1 5],'InitialMagnification','fit'); colormap(jet); colorbar
(F2=immagine da trasformare, [-1 5]=range di tonalità da mostrare, InitialMgnification e
Fit=riduce la grandezza originale allo schermo di visualizzazione).
29
La visualizzazione risulta diversa dalla forma che ci si aspettava: il campionamento è
grossolano e la continua si trova agli angoli dell’immagine.
3)Miglioramento visualizzazione
3.1) Raffinare il campionamento, tramite lo zero-padding fino a rendere l’immagine di
256x256
F = fft2(f,256,256);
imshow(log(abs(F)),[-1 5]); colormap(jet); colorbar
3.2)Portare la componente continua al centro dell’immagine, usando la funzione fftshift ,
che sposta i quadranti della trasformata.
F = fft2(f,256,256);F2 = fftshift(F);
imshow(log(abs(F2)),[-1 5]); colormap(jet); colorbar
Una applicazione della FFT implementabile con MATLAB è quella del calcolo della
trasformata della risposta all’impulso di un filtro. La funzione freqz2 calcola e mostra la
risposta in frequenza di un filtro, il cui tipo va specificato attraverso la funzione fspecial. Ad
esempio se si vuole visualizzare la risposta in frequenza di un filtro gaussiano che ha
comportamento passa-basso, si dovrà specificare
h = fspecial('gaussian');
freqz2(h)
Si è già parlato dell’importanza della FFT nel calcolo di convoluzioni di immagini per
grandi operatori, sfruttando la proprietà di Fourier che il prodotto di risposte in frequenza
corrisponde alla convoluzione nel dominio delle spazio. L’esempio di seguito esegue
proprio questo tipo di operazione tra due matrici A e B rispettivamente MxN e PxQ
A = magic(3);
%creare le matrici
B = ones(3);
A(8,8) = 0;
B(8,8) = 0;
C = ifft2(fft2(A).*fft2(B));
C = C(1:5,1:5);
C = real(C)
% rendere le loro dimensioni una potenza
di 2 attraverso zero-padding
%estrazione della porzione di immagine pulita dagli zeri
%visualizzazione parte reale
30
Altra possibilità di impiego della FFT si è detto essere la pattern-recognition . Le righe di
codice seguenti mostrano l’implemetazione MATLAB di questo processo, già descritto in
precedenza.
bw = imread('text.png');
%estrae da ‘bw’ l’immagine della lettera “a”
a = bw(32:45,88:98);
imshow(bw);
%mostra l’immagine aprendo una nuova
finestra
figure, imshow(a);
C = real(ifft2(fft2(bw) .* fft2(rot90(a,2),256,256)));
%effettuo rotazione del target, poiché la
correlazione e la convoluzione differiscono
nel fatto chela correlazione necessita che il
target sia ruotato di 180°
figure, imshow(C,[])
%mostra imagine dopo
adattandola al display
averla
scalata
max(C(:))
ans =
%trova il valore del pixel di intensità maggiore
per individuare l’elemento nell’immagine
68
%fissa una soglia minore del max
%mostra solo i pixel sopra la soglia per una
individuazione più semplice degli elementi
trovati
thresh = 60;
figure, imshow(C > thresh)
4.5.Altri tipi di trasformata
4.5.1. Trasformata Coseno Discreta
La Trasformata Coseno Discreta, o DCT, fornisce un metodo di rappresentazione di una
immagine come somma di sinusoidi di ampiezza e frequenza diversa.
La DCT unidimensionale è definita dalla seguente formula:
F (u ) =
N −1
2
( 2 x + 1)uπ
C (u )∑ f ( x ) cos[
]
N
2N
x =0
31
a DCT bidimensionale, per una matrice quadrata di NxN elementi, si ottiene applicando la
definizione precedente a ogni riga e successivamente a ogni colonna della matrice. Per ogni
riga e colonna vale la formula precedente, quindi risulta:
N−1 N−1
2
(2x +1)uπ
(2y +1)vπ
F(u, v) = C(u)C(v)∑∑ f ( y, x) cos[
]cos[
]
N
N
N
2
2
x=0 y=0
Anche per questo tipo di trasformata, il coefficiente F(0,0), che rappresenta la frequenza
spaziale 0 in orizzontale e verticale, è chiamato coefficiente “DC”. Tutti gli altri coefficienti
sono “AC”.
Una caratteristica di questa trasformata è che le informazioni circa l’immagine sono
concentrate in pochi coefficienti relativamente scorrelati, rilevati come variazioni di
informazione tra un'area e quella contigua trascurando le ripetizioni; è per questo che questo
tipo di trasformata fornisce la base per il semplice algoritmo di compressione con perdita
delle immagini, il JPEG. Altra caratteristica della DCT è quella di offrire una
rappresentazione dei dati nel dominio delle frequenze spaziali, permettendo di ridurre la
precisione dei coefficienti all'aumentare della frequenza, sfruttando la sensibilità dell'occhio
umano: paragonata con la DFT, si osserva, quindi, che lo spettro della DFT è più diffuso
dello spettro della DCT e che la DCT concentra le informazioni nelle basse frequenze.
La DTC è una trasformata simile alla trasformata discreta di Fourier (DFT), ma fa uso solo
di numeri reali. È equivalente a una DFT di lunghezza circa doppia, che opera su funzioni
reali e pari (dato che la trasformata di Fourier di una funzione reale e pari è reale e pari a sua
volta), dove in alcune varianti l'input e/o l'output sono traslati di mezzo campione.
La DFT e' uno sviluppo del segnale su una base di coseni e seni, mentre la DCT e' uno
sviluppo su una base di soli coseni. E' possibile sviluppare un segnale reale in soli coseni se
e solo se la sua trasformata di Fourier è reale pura, cioè se il segnale è pari (simmetrico
32
attorno a 0). E' possibile sviluppare un segnale qualunque in serie di soli coseni se il
dominio del segnale è raddoppiato mediante una “specchiatura” attorno a 0 dello stesso
segnale:
La DCT è una trasformata di tipo invertibile e l’inversa di una matrice MxN è interpretabile
come la somma di MN funzioni chiamate funzioni base;per una immagine 8x8 le 64 basi
sono rappresentate nella figura sotto dove è possibile osservare come le frequenze
orizzontali migliorano da sinistra e destra e le verticali dall’altro al basso.
Ci sono due modi per operare con questa trasformata in MATLAB, il primo metodo implica
l’utilizzo della funzione dct2 che utilizza un algoritmo basato sull’impiego della FFT, ed è
utilizzata nel caso in cui voglia trasformare elementi grandi. Il secondo metodo consiste
nell’utilizzare la matrice di trasformazione attraverso l’utilizzo della funzione dctmtx, ed è
maggiormente impiegato per trasformare piccoli elementi, come matrici 8x8 o 16x16.
33
4.5.2.Trasformata di Radon
La trasformata di Radon (RT) è una trasformata integrale la cui inversa (detta
antitrasformata di Radon) è utilizzata per ricostruire immagini bidimensionali a partire dai
dati raccolti nel processo di diagnostica medica di tomografia assiale computerizzata (TAC).
Conoscere la trasformata di Radon di un oggetto permette di ricostruirne la struttura: il
"teorema della proiezione" infatti assicura che se abbiamo un numero infinito di proiezioni
monodimensionali di un oggetto fatte da un numero infinito di angoli diversi, possiamo
ricostruire perfettamente la geometria dell'oggetto originale e il processo di ricostruzione
consiste appunto nel calcolare l'antitrasformata di Radon.
La trasformata di Radon R della funzione f(x,y) è definita come:
R(m, q)[ f ( x, y)] = ∫ f ( x, mx + q)dx = ∫∫ f ( x, y)δ [ y − (mx + q)]dxdy = U (m, q)
dove m è la pendenza (o coefficiente angolare) di una retta e q è la sua quota (o intercetta).
In MATLAB questo tipo di trasformata è affidata alla funzione radon, che calcola appunto
le proiezioni dell’intensità di una immagine lungo linee radiali orientate con un angolo
specificato, dove per proiezioni si intendono integrali di linea. Questa funzione calcola gli
integrali da diverse sorgenti lungo degli assi o linee parallele distanti tra loro 1 pixel. Per
rappresentare una immagine, quindi, questa funzione acquisisce proiezioni dell’immagine
da diverse angolazioni.
Ad esempio, l’integrale di linea di f(x,y) nella direzione verticale altro non è che la
proiezione lungo l’asse x e l’integrale di linea nella direzione orizzontale quella lungo l’asse
y. Le proiezioni dell’immagine precedente risultano essere
34
Di solito la trasformata di Radon di un f(x,y) coincide con la proiezione lungo l’asse y.
La funzione radon, applicata ad una immagine I e richiesta per un angolo theta, viene
utilizzata con la seguente sintassi
[R,xp] = radon(I,theta);
e restituirà un vettore colonna R, contenente le trasformate per ogni angolo specificato.
La rotazione è calcolata sulla base del valore centrale dell’immagine, che si determina
attraverso una linea di codice del tipo
floor((size(I)+1)/2);
Le linee di codice seguenti creano una immagine quadrata e operano la RT a 0° e 45°.
I = zeros(100,100);
I(25:75, 25:75) = 1;
imshow(I)
[R,xp] = radon(I,[0 45]);
figure; plot(xp,R(:,1)); title('R_{0^o} (x\prime)')
35
figure; plot(xp,R(:,2)); title('R_{45^o} (x\prime)')
Nel caso in cui si voglia visualizzare la trasformata lungo un intervallo di angoli, ad
esempio tra 0° e 180° con incrementi di 1°, si procede definendo theta = 0:180 e in questo
caso è possibile visulizzare la trasformata in forma di immagine, come descritto nelle righe
di codice seguenti
theta = 0:180;
[R,xp] = radon(I,theta);
imagesc(theta,xp,R);
title('R_{\theta} (X\prime)');
xlabel('\theta (degrees)');
ylabel('X\prime');
set(gca,'XTick',0:20:180);
colormap(hot);
colorbar
Una applicazione di questa trasformata è il riconoscimento di linee rette, magari derivanti
da disturbi di acquisizione, all’interno di una immagine, come effettuato dalla trasformata
di Haugh. Per compiere questa operazione si procede, dopo aver trasformato l’immagine in
immagine binaria, ad effettuare la trasformata e a rilevare nell’immagine derivante i picchi
di intensità che individuano linee rette.
36
Per quanto riguarda l’antitrasformata di Radon, la sua espressione matematica è data dalla
seguente formula
[ f ( x, y)] =
1 d
H[U (m, y − mx)]dm
2π ∫ dy
dove H è la trasformata di Hilbert.
Questa operazione è utilizzata per ricostruire immagini precedentemente trasformate,
attraverso la funzione iradon secondo la formulazione IR = iradon(R,theta);dove la R indica
la trasformata di una immagine I, e quindi l’insieme delle proiezioni di quell’immagine.
Nella maggior parte dei campi di applicazione, però non si dispone di immagini da cui
estrarre le proiezioni, come per esempio nel caso delle applicazioni alla tomografia, nelle
quali le proiezioni sono formate misurando l'attenuazione della radiazione che passa
attraverso uno spazio da diverse angolazioni.
L'immagine originale può essere vista come una sezione trasversale di questo spazio, in cui
l'intensità valori rappresentano la densità del campione. Le proiezioni sono raccolte da un
apparato specifico e da queste viene ricostruita l’immagine del campione attraverso la
funzione iradon.
5.FILTRAGGIO
5.1.DOMINIO DELLO SPAZIO
Il filtraggio delle immagini è una operazione “spaziale” nel senso che considera ogni
singolo pixel dell’immagine in ingresso come elemento di una zona di pixel con
caratteristiche simili e che quindi fornisce una immagine in uscita in cui ogni pixel assume
un valore che è stato dedotto attraverso l’applicazione all’immagine di ingresso di algoritmi
37
che attuano una analisi dei pixel circostanti: la trasformazione dei pixel avviene
manipolando direttamente i pixel stessi, mediante l’uso di maschere (o filtri spaziali), dove
la maschera è una matrice di pixel nell’intorno del pixel considerato.
Il funzionamento del filtro spaziale è quello di muovere la finestra sull'immagine di un pixel
sulle colonne e poi sulle righe e per ogni posizione applicare una determinata relazione
funzionale ai valori dei pixel ricadenti la sottoimmagine.
Esistono due tipi di filtri spaziali: i filtri lineari o di convoluzione e i filtri non lineari. Per i
filtri lineari, la risposta dei filtri spaziali è data da una combinazione lineare dei pixel della
sottoimmagine, utilizzando opportuni coefficienti che definiscono il filtro: si compie una
operazione di convoluzione che, nell’immagine seguente è rappresentata per una maschera
3x3.
Esempio di filtro lineare è il filtro di smoothing (passa-basso) che riduce il rumore presente
nelle immagini, ma tende a peggiorare la definizione dei dettagli, in particolare dei bordi.
Il filtro di smoothing tende ad eliminare i dettagli irrilevanti le cui dimensioni siano
tipicamente più piccole della maschera, attenuando ad esempio gli effetti dell’aliasing,
oppure come riduzione del rumore. Quest'ultima applicazione è la più investigata, ma ha lo
svantaggio di ridurre oltre al rumore anche le variazioni rapide di livello di grigio associate
ai contorni degli oggetti (edge) e non al rumore. E’ possibile compiere questo tipo di
filtraggio o tramite Box Filter, che ha proprietà di preservare la media, intesa come il valore
di luminanza medio nello spazio nel caso di segnali bidimensionali, o anche attraverso una
maschera weighted average ossia a media pesata, che avrà valori che decrescono lentamente
man mano che ci si allontana dal centro della maschera e in modo proporzionale alla
distanza euclidea.
Ricapitolando, le operazioni che si compiono al fine di filtrare linearmente una immagine
sono operazioni di convoluzione (ogni pixel in out è la somma pesata di quelli circostanti il
corrispondente pixel in ingresso attraverso una matrice (ruotata di 180 gradi)) e di
correlazione (la sola differenza con la precedente è che la matrice dei pesi non è ruotata).
38
Queste operazioni sono compiute attraverso la funzione imfilter. Ad esempio per filtrare
linearmente tramite una operazione di media si possono usare le seguenti linee di codice che
descrivono il filtraggio dell’ immagine “coins.png” attraverso un filtro 5 per 5 di pesi
uguali.
I = imread('coins.png');
h = ones(5,5) / 25;
I2 = imfilter(I,h);
imshow(I), title('Original Image');
figure, imshow(I2), title('Filtered Image')
E’ inoltre possibile utilizzare filtri diversi da quelli lineari tramite la funzione fspecial che
produce alcuni tipi di filtri predefiniti, che saranno poi applicabili all’immagine tramite la
funzione imfilter. Nell’esempio viene utilizzata la maschera ‘unsharp’ che esalta i contorni e
rende l’immagine più nitida.
I = imread('moon.tif');
h = fspecial('unsharp');
I2 = imfilter(I,h);
imshow(I), title('Original Image')
figure, imshow(I2), title('Filtered Image')
I filtri non lineari considerano sempre una matrice di pixel per operare la trasformazione,
ma non utilizzano coefficienti: mediante metodi statistici o formule matematiche ricavano
un nuovo valore per l'intensità del pixel considerato, a partire dall'insieme dei valori della
matrice.
Esempio di filtro non lineare è il filtro mediano che opera sostituendo il valore del pixel
considerato con la mediana dei valori della matrice, riducendo così il rumore nelle
immagini, ma preservando maggiormente i dettagli e i contorni del filtro di sfocatura. In
particolare viene effettuato un ordinamento dei pixel all'interno di una regione e poi viene
prelevato il valore centrale e assegnato al pixel da processare. L'istruzione di Matlab che
implementa questo filtro è medfilt2. Questa istruzione può risultare dannosa nel caso di
rumore additivo gaussiano o uniforme, è invece molto efficace nel caso di rumore
impulsivo, aggiunto da Matlab con l'opzione salt pepper in imnoise. In tal modo infatti, se il
rumore sale e pepe si manifesta su un numero di pixel non elevatissimo, i punti neri (pepe) e
quelli bianchi (sale) appariranno come delle code (outliers) nelle sottoimmagini e saranno
sostituiti dal valore mediano che tipicamente sarà determinato dai pixel non rumorosi
purché in numero elevato (questo filtraggio verrà analizzato meglio nel paragrafo dedicato).
Un altro esempio di filtro non lineare è il filtro di sharpening che evidenzia le differenze
d'intensità tra i pixel, in particolare i contorni degli oggetti, riducendo gli altri dettagli a
sfondo, ma che, d’altra parte, enfatizza anche il rumore nell'immagine, che va eliminato
prima dell’applicazione del filtro. Se per arrotondare i dettagli (smoothing) ci si è serviti
39
della media che è dopotutto una forma di integrazione, allora per esaltarli (sharpening) ci si
servirà di alcune forme di derivazione. Nei filtri derivativi l’intensità della risposta è
proporzionale al grado di discontinuità dell’immagine nel punto in cui è applicata. La
differenziazione accentua i bordi e altre discontinuità (come il rumore) e trascura le aree a
lenta variazione dei livelli di grigio. Possono essere considerati filtri derivativi del primo e
del secondo ordine e per ciascuno di essi sono possibili molteplici approssimazioni discrete
sotto forma di differenze finite. I filtri derivativi del secondo ordine sono più sensibili alle
brusche transizioni e sono utilizzati per l’enhancement dei dettagli (come i filtri Laplaciani).
Quelli del primo ordine sono indicati nei processi di edgedetection e extraction ( filtri di
contorno orizzontale e verticale).
Ulteriori esempi di filtri non lineari sono forniti dai filtri di dilatazione ed erosione nei
quali il valore del pixel viene sostituito rispettivamente con il valore massimo e minimo
della matrice (il primo espande gli oggetti presenti nell'immagine, il secondo li riduce) e dai
filtri di differenziazione (Roberts, Sobel, Fase), basati sulla determinazione del gradiente
dell'intensità, e che evidenziano i contorni.
5.2.DOMINIO DELLE FREQUENZE
Il nucleo del filtraggio in frequenza consiste nel modificare in qualche modo la trasformata
dell’immagine e nel prendere poi la sua antitrasformata come immagine elaborata in uscita;
si procede a questo cambio di dominio grazie agli ampi vantaggi e della semplicità
dell’elaborazione in frequenza che sono stati presentati nel paragrafo dedicato.
Ovviamente anche in frequenza è possibile compiere operazioni di smoothing e sharpening:
per quanto concerne lo smoothing, che serve ad eliminare le brusche transizioni ed il
rumore; questo corrisponde nel dominio della frequenza ad attenuare le componenti ad alta
frequenza. Esistono molti tipi di lowpass filter: i principali sono quello ideale, quello di
Butterworth e quello Gaussiano. Per quanto riguarda il filtraggio ideale, tutte le frequenze
all’interno di un cerchio di raggio D0 sono trasferite senza attenuazione, tutte quelle esterne
sono annullate; l’idealità comporta impossibilità di una realizzazione fisica (sebbene
possano realizzarsi in digitale con l’ausilio di un computer). Al variare della frequenza di
taglio D0 il filtro permetterà il passaggio di diverse porzioni dello spettro di potenza
dell’immagine, anche se in ogni caso questo tipo di filtraggio provoca un forte fenomeno di
sfocatura ad anello. L’origine dell’effetto di ringing è legato alla funzione antitrasformata
del filtro lowpass ideale.
40
Il filtraggio alla Butterworth, invece, per basso ordine, non determina un taglio netto in
frequenza, in quanto la sua funzione di trasferimento non esibisce una discontinuità brusca
come mostrato in figura e quindi non presenta l’effetto di ringing.
Altro filtraggio è quello Gaussiano che ha la caratteristica che la sua antitrasfomata è ancora
una gaussiana. Le prestazioni di questo tipo di filtro sono simili a quelli alla Butterworth per
basso ordine, ma di più facile realizzazione.
Per quanto riguarda invece lo sharpening, filtraggio passa-alto per correggere le brusche
variazioni, si utilizzano lo stesso tipo di filtri con caratteristiche e problematiche molto
simili ai corrispondenti low-pass.
41
In MATLAB filtraggio in frequenza consiste nel creare un filtro bidimensionale a partire da
uno monodimensionale creato attraverso le funzioni del Signal Processing Toolbox (la
bidimensionalità è richiesta dalla bidimensionalità dell’immagine). In generale, ogni filtro
che può essere utilizzato nel toolbox di matlab è un filtro FIR. Questo tipo di filtri
forniscono la base ideale per l’elaborazione matlab delle immagini, in quanto sono
facilmente rappresentabili come matrici di coefficienti, esistono già metodi certi di progetto,
sono facili da implementare ed hanno fase lineare, quindi prevengono le distorsioni. Si
chiamano filtri Finite Impulse Response (FIR) tutti i filtri caratterizzati da una risposta
all'impulso con un numero finito di campioni.
Essi si realizzano mediante l'operazione di convoluzione: per ogni campione di uscita
prodotto viene calcolata una somma pesata di un numero finito di campioni dell'ingresso.
Primo tipo di filtraggio è quello della trasformazione di frequenza. Le trasformazioni in
frequenza consistono nel trasformare un filtro FIR ad una dimensione in uno a due
dimensioni e conservano caratteristiche come la banda di transizione e il ripple. Questo
metodo utilizza una matrice di trasformazione e la funzione ftrans2, che crea filtri a
simmetria centrale. Questo metodo fornisce molto spesso buoni risultati, perché permette
partire da un filtro ad una dimensione, che è sicuramente più facile da progettare, per
arrivare ad un filtro a due dimensioni, necessario per processare una immagine,
mantenendone la forma e le caratteristiche richieste, come mostra la simulazione.
42
Altro tipo di filtraggio è quello basato sul campionamento in frequenza. Questo metodo crea
un filtro basato su una risposta in frequenza desiderata. Data una matrice di punti che
definiscono la forma della risposta in frequenza, questo metodo crea un filtro la cui risposta
in frequenza passa attraverso questi punti. Il compito è svolto dalla funzione fsamp2, ed è un
metodo che opera direttamente su un filtro a due dimensioni.
Come si può osservare, lo svantaggio di questo metodo è la presenza di un forte ripple in
banda derivante da discontinuità nette nel filtro ideale.
Un filtraggio migliore rispetto a quello appena descritto è quello che si ottiene tramite
finestramento: si moltiplica la risposta all’impulso ideale per una funzione a finestra in
43
modo da ottenere il filtro corrispondente. In MATLAB sono presenti due funzioni attraverso
le quali è possibile utilizzare il metodo del finestramento, fwind1 e fwind2. (La funzione
fwind1 crea un filtro bidimensionale usando una finestra bidimensionale derivante da una o
due finestre monodimensionali precedentemente specificate, invece la funzione fwind2 crea
un filtro bidimensionale utilizzando direttamente una finestra bidimensionale
precedentemente specificata).
Le funzioni introdotte per gli ultimi due metodi di progetto di filtro in frequenza richiedono
di conoscere a priori la risposta in frequenza desiderata; è possibile creare la matrice
corrispondente alla risposta in frequenza desiderata attraverso la funzione freqspace.
L’immagine in basso mostra il filtro 11 per 11 creato attraverso la funzione fwind1 dalla
funzione desiderata Hd (la funzione hamming ne crea la finestra monodimensionale
corrispondente).
Hd = zeros(11,11); Hd(4:8,4:8) = 1;
[f1,f2] = freqspace(11,'meshgrid');
mesh(f1,f2,Hd), axis([-1 1 -1 1 0 1.2]), colormap(jet(64))
h = fwind1(Hd,hamming(11));
figure, freqz2(h,[32 32]), axis([-1 1 -1 1 0 1.2])
6.OPERAZIONI MORFOLOGICHE
6.1.Dilatazione ed erosione_ Questo tipo di operazioni applicano ad una immagine
in ingresso un dato strutturato (matrice binaria) in modo da creare in uscita una immagine
della stessa dimensione di quel dato. Per elaborare una immagine spesso si utilizzano
processi di scheletrizzazione (ridurre tutti gli elementi di una immagine ad una linea, usando
la funzione bwmorph) e determinazione dei bordi (determinare i pixel dei perimetri degli
oggetti dell’immagine tramite la funzione bwperim). E’ possibile dilatare ripetutamente una
immagine in modo tale da rendere i suoi contorni uguali a quelli di una immagine maschera:
si sta effettuando una ricostruzione morfologica. E’ anche possibile effettuare una
44
trasformazione di distanza: misurare la distanza tra due punti di una immagine. (La funzione
bwdist misura la distanza di un pixel 0 e il pixel 1 più vicino in una immagine binaria.
6.2.Individuare e misurare oggetti all’interno di una immagine_ Un
componente connesso in una immagine binaria è un insieme di pixel che forma un gruppo
connesso (vedi “pixel connettivity”).
L’individuazione di oggetti all’interno dell’immagine è il processo di identificazione di un
gruppo connesso di pixel e l’assegnazione ad ogni gruppo di un identificativo, come
mostrato in figura.
La funzione bwconncomp rileva il componente connesso. Ad esempio
cc = bwconncomp(BW)
cc =
Connectivity: 8
ImageSize: [8 9]
NumObjects: 3
PixelIdxList: {[6x1 double] [6x1 double] [5x1 double]}
(La voce PixelIdxList identifica I pixel dei singoli componenti connessi)
45
Per visualizzare i componenti connessi si usa la funzione labelmatrix attraverso la quale ad
ogni oggetto viene associato un colore nella matrice colormap associata.
labeled = labelmatrix(cc);
Si usa la funzione label2rgb per scegliere i parametri attraverso i quali definire la colormap.
RGB_label = label2rgb(labeled, @copper, 'c', 'shuffle');
imshow(RGB_label,'InitialMagnification','fit')
6.3.Selezionare un singolo oggetto in una immagine binaria_ E’ possibile
utilizzare la funzione bwselect per selezionare un singolo oggetto, specificando dei pixel
nell’immagine di ingresso, selezionati anche da mouse; la funzione restituirà una immagine
binaria che include solo quegli oggetti dell’immagine di ingresso che contengono i pixel
specificati.
6.4.Operazioni tramite uso di Look up Table_
Alcune operazioni su
immagini sono implementabili più facilmente tramite l’uso di look up table: vettori colonna
i cui elementi sono tutti i possibili valori restituibili come risultato di combinazione dei
pixel nell’area che si sta analizzando. Si utilizzano aree di 2per2 o 3per3 e la funzione che
crea look up table su aree di queste dimensioni è la makelut. Una volta creata, è possibile
utilizzare la LUT attraverso la funzione applylut . L’esempio di seguito crea una funzione
che restituisce 1 se tre o più pixel in un’area 3per3 è 1, altrimenti restituisce 0 e la passa in
ingresso alla makelut, insieme all’indicazione della grandezza dell’area da considerare;
successivamente compie l’operazione voluta tramite applylut.
f = @(x) sum(x(:)) >= 3;
lut = makelut(f,3);
BW1 = imread('text.png');
BW2 = applylut(BW1,lut);
imshow(BW1)
46
figure, imshow(BW2)
7.ANALIZZARE LE IMMAGINI
7.1.Acquisire informazioni sul valore dei pixel e sulla statistica delle
immagini_ Per conoscere il valore di uno o più pixel individuati tramite coordinate o da
mouse viene utilizzata la funzione impixel .
Per creare un profilo di intensità dell’intera immagine si usa improfile, o per creare un
istogramma della distribuzione dell’intensità si usa imhist.
7.2.Eliminare il rumore da una immagine_ Il rumore in una immagine non è
altro che una variazione nel valore dei pixel che non rende l’immagine fedele alla realtà; ci
sono diversi modi attraverso i quali il rumore è introdotto in dipendenza da come
l’immagine è creata, ed è possibile simulare gli effetti di questi tipi di rumore aggiungendo
ad una immagine del rumore prodotto attraverso la funzione imnoise.
L’eliminazione del rumore può avvenire con l’impiego di un filtraggio di tipo lineare ad
opera di filtri di media o gaussiani. Un filtro di media è molto impiegato nel caso di rumore
a grani: ogni pixel dell’immagine filtrata deriva dalla media dei pixel circostanti il
corrispondente pixel in ingresso, quindi ogni variazione locale data da rumore a grani è
ridotta. Una funzione utilizzata per il filtraggio di media e la filter2.
K = filter2(fspecial('average',3),J)/255
Altro tipo di filtraggio che si può impiegare è quello mediano, che è simile a quello di
media, ma differisce per il fatto che ogni pixel deriva si dalla media di quelli circostanti, ma
da una media di tipo statistico, che risulta quindi essere molto meno sensibile a valori rari, e
quindi questo tipo di filtro è migliore nell’eliminare questi valori poco frequenti
preservando la nitidezza dell’immagine. Questo tipo di filtri è implementato dalla funzione
medfilt2.
L = medfilt2(J,[3 3])
E’ possibile inoltre utilizzare un filtraggio adattativo: la funzione wiener2 applica un filtro
Wiener (tipo di filtro lineare) adattandolo all’immagine alla quale lo applica, in particolare
adattandolo alla varianza locale dell’immagine: se la varianza è grande, il filtro attenua
poco, altrimenti tanto. Questo metodo risulta essere più selettivo rispetto quello lineare, e
non sporca i contorni e altre componenti di alta frequenza dell’immagine. Questo tipo di
filtraggio è usato soprattutto nel caso di rumore bianco (tipo rumore gaussiano).
K = wiener2(J,[5 5])
47
8.UTILIZZARE UNA REGIONE DI INTERESSE
(ROI)
8.1.ROI_ La Regione di Interesse è una porzione di immagine sulla quale si vuole
operare che si definisce attraverso l’utilizzo di una maschera binaria di dimensioni pari alla
zona da evidenziare. Si possono individuare più zone di interesse sulla base di principi
geometrici, o anche di intensità. Data una immagine, per creare la maschera binaria,
composta di 1 all’interno del ROI e 0 all’esterno,si utilizza il metodo createMask, e quindi
si usa una delle funzioni seguenti sulla base del criterio di selezione per l’estrazione del
ROI: impoint, imline, imrect, imellipse, impoly o imfreehand. Ad esempio, supponendo di
voler rilevare una porzione di immagine circolare si procederà così
img = imread('pout.tif');
h_im = imshow(img);
e = imellipse(gca,[55 10 120 120]);
BW = createMask(e,h_im);
Si può anche creare un ROI senza una immagine di riferimento, tramite la funzione
poly2mask, specificando i vertici del ROI in due vettori e la dimensione della maschera che
si vuole ottenere. Ad esempio,le linee di codice seguenti creano una maschera binaria a
prescindere dall’immagine alla quale applicarla.
c = [123 123 170 170];
r = [160 210 210 160];
m = 291; % altezza immagine in uscita
n = 240; % larghezza immagine in uscita
BW = poly2mask(c,r,m,n);
figure, imshow(BW)
Altro metodo per la creazione di ROI è quello basato sull’intensità di colore: si utilizza la
funzione roicolor per creare un ROI sulla base di un range di intensità di colore.
8.2.Filtrare una ROI_ Filtrare una ROI, non è altro che filtrare una zona di
immagine definita dall’applicazione di una maschera, utilizzando la funzione roifilt2
specificando la scala di grigi in ingresso da filtrare, la maschera che definisce il ROI e il tipo
di filtro.
48