Metodi Bayesiani per problemi inversi e applicazioni in Matlab all'elaborazione delle immagini

1 Perturbazioni delle immagini
digitali
Le immagini digitali sono funzioni di due o tre variabili (rispettivamente nel
caso di immagini 2D o 3D) prodotte da sistemi sici di registrazione dei segnali,
degli strumenti in grado trasformare le radiazioni rilevate (fotoni, raggi, onde, ecc)
in radiazioni contenenti informazioni, e di quanticarne l’intensit a.
Le immagini cos  ottenute sono composte da pixel (da picture element) nel caso
2D, o da voxel (da volume element) nel caso 3D, che contengono le informazioni
riguardo un oggetto sico; ad ognuno  e assegnato un colore che rappresenta il
livello di espressione della quantit a sica misurata (o di una sua funzione), in
corrispondenza di quella posizione spaziale secondo una determinata gradazione.
Formazione delle immagini e perturbazioni
Indicato conX il dominio spaziale nel quale si deve risolvere il problema, che
sar a in generale un sottoinsieme di R
2
o di R
3
, e si suppone di aver ssato un
opportuno sistema di riferimento cartesiano inX, in modo tale da poter rappre-
sentare ogni punto (cio e ogni posizione) con un vettore x. Si assume inoltre che i
vettori siano vettori colonna nito-dimensionali.
I processi di formazione e registrazione delle immagini sono descritti mediante
trasformazioni  fra due spazi funzionali su due dominiX eY
 :F(X)!G (Y ) (1.1)
1
1 Perturbazioni delle immagini digitali
conF(X) spazio dell’oggetto eG(Y ) spazio dell’immagine. Gli elementi diF(X)
sono funzioni suX che rappresentano la distribuzione ‘vera’ della propriet a sica
che si vuole indagare nel dominio X (incognita del problema). Gli elementi di
G(Y ) sono funzioni suY che rappresentano invece la distribuzione della quantit a
eettivamente misurata nell’osservazione dell’oggetto (dati del problema). Il do-
minioY  e anch’esso un sottoinsieme di R
n
e pu o o meno coincidere con X.
Ad esempio, nel processo di formazione delle immagini radiograche di un corpo, il
dominioX  e un volume tridimensionale contenente il corpo mentre il dominio Y  e
una parte di piano. In tal caso la funzione incognitaf(x)2F(X) rappresenta la
distribuzione (incognita) della densit a del corpo mentre g(y)2G(Y ) rappresenta
la distribuzione misurata (e quindi nota) dell’attenuazione della radiazione che ha
attraversato il corpo (cio e il dato).
Le immagini che si ottengono tramite questo processo non sono per o fedeli alla
realt a, per due cause principali:
 la presenza di disturbi nel processo di formazione dell’immagine, che pos-
sono essere dovuti a diverse cause, come per esempio il moto relativo tra
la macchina fotograca e i soggetti della foto (che causa immagini mosse
e sfuocature), le aberrazioni ottiche (soprattutto nel caso di microscopi e
telescopi) o le turbolenze atmosferiche o polveri; si parla in questo caso di
blurring;
 la comparsa di disturbi nel processo di registrazione dell’immagine dovuti
per esempio alla discretizzazione o ad errori di musura e di conteggio; si
parla in questo caso di noise.
In genere il blurring  e dato da un processo deterministico e sono disponibili model-
li matematici per descriverlo e risalire all’immagine originale. Nel caso del noise
invece, si tratta di un processo stocastico che non  e facile determinare.
2
1.1 Blurring
1.1 Blurring
Nella formazione delle immagini ottiche, il segnale in ingresso  e una distribuzione
sconosciuta di radiazioni luminose, chiamate oggetto e denotate con f
(0)
(x), dove
x  e un vettore bidimensionale, cio e del tipo x = (x
1
;x
2
)
T
, oppure tridimensionale
cio e del tipo x = (x
1
;x
2
;x
3
)
T
, che produce una distribuzione radiale nel dominio
del sistema ottico usato per registrarla. Questa si chiama immagine priva di noise
di f
(0)
(x) e si denota con g
(0)
(x), che  e quindi il segnale in uscita.
Nel caso in cui il processo di formazione dell’immagine si possa modellizzare, come
spesso accade, in modo lineare, si ha che
g
(0)
(x) =
Z
K(x;x
0
)f
(0)
(x
0
)dx
0
(1.2)
doveK(x;x
0
)  e la funzione di dispersione puntuale ( Point Spread Function, PSF)
del sistema lineare, anche chiamata funzione di allargamento dei punti del sistema
oppure risposta in impulso del sistema. A volte si usa anche il termine funzio-
ne di diusione dei punti per indicare K(x;x
0
), perch e rappresenta l’immagine
prodotta da una sorgente di luce puntiforme localizzata nel puntox
0
, solitamente
modellizzata mediante una distribuzione (x) di Dirac concentrata inx
0
. Infatti se
la sorgente del punto  e data da f
(0)
(x
00
) = (x
00
 x
0
) si ottieneg
(0)
(x) =K(x;x
0
).
La PSF d a quindi l’immagine di sorgenti puntiformi di intensit a unitaria e fornisce
un’immagine in cui i dettagli sono meno nitidi rispetto all’immagine originaria;
per questo l’eetto di tale funzione  e chiamato blurring (annebbiamento) e l’im-
magine g
(0)
(x)  e chiamata versione con blurring di f
(0)
(x) (si pensi ad esempio
alla deformazione introdotta da un occhio miope guardando una stella, che la fa
vedere come una sorta di macchia).
In generale, la PSF si pu o ottenere risolvendo il problema diretto associato al pro-
cesso di formazione dell’immagine. Nel caso di un sistema ottico questo consiste
nel calcolare la propagazione della luce da un punto sorgente nel dominio ad un
3
1 Perturbazioni delle immagini digitali
punto del codominio attraverso gli elementi del sistema (lenti, specchi, ecc). Se
non si dispone di una rappresentazione analitica, neanche approssimata, del siste-
ma diretto, la PSF pu o essere misurata generando un punto sorgente, muovendolo
nel dominio e valutando le immagini prodotte nelle varie posizioni.
Molti sistemi di formazione di immagini sono (o vengono approssimati come) iso-
planatici, ossia invarianti per traslazioni; ne segue che anche la loro PSF  e inva-
riante rispetto alle transazioni, cio e la K(x;x
0
) di una sorgente che si trova inx
0
si pu o anche vedere come la traslata in x
0
della PSF di una sorgente che si trova
nell’origine, cio e K(x;x
0
) =K(x x
0
; 0). Ne segue cheK(x;x
0
)  e funzione della
dierenza x x
0
e si pu o quindi scrivere K(x;x
0
) =K(x x
0
). Questo implica che
la forma dell’immagine di un punto luminoso non dipende dalla sua posizione: la
PSF, che dunque si pu o scrivere semplicemente come K(x), rappresenta l’immagi-
ne di una sorgente puntiforme di intensit a unitaria posta al centro dell’immagine.
Una PSF spazio-invariante pu o quindi essere determinata guardando l’immagine
di un singolo punto sorgente. In questo caso particolare si ha inoltre che la (1.2)
diventa
g
(0)
(x) =
Z
K(x x
0
)f
(0)
(x
0
)dx
0
(1.3)
o analogamente, per denizione di convoluzione,
g
(0)
=K f
(0)
: (1.4)
La funzioneK(x)  e sempre chiamata PSF e la sua trasformata di Fourier
^
K(!)  e
chiamata funzione di trasferimento (TF) del sistema di immagini.
In termini di trasformata la (1.4) diventa quindi
^ g
(0)
(!) =
^
K(!)
^
f
(0)
(!): (1.5)
Queste due funzioni giocano un ruolo molto importante. La PSF fornisce la rispo-
sta del sistema per ogni punto sorgente, dovunque tale punto si trovi all’interno
4
1.1 Blurring
del dominio. La TF dice quanto segnale di una frequenza ssata  e propagato
attraverso il sistema lineare; il blurring pu o essere visto allora come una sorta di
ltrazione della frequenza.
La descrizione di un sistema di formazione di immagini tramite la sua PSF viene
anche detta descrizione nel dominio spaziale, mentre la sua descrizione mediante
la TF viene detta descrizione nel dominio delle frequenze. Spesso  e utile conside-
rare il sistema da questo secondo punto di vista, ad esempio quando la PSF  e a
banda limitata, cio e nulla al di fuori di un dominio.
In generale si dice che un sistema per immagini  e limitato se la sua PSF  e limita-
ta. La banda B della PSF  e anche chiamata banda del sistema. Negli esempi che
seguiranno il sistema sar a limitato oppure approssimatamente limitato.
Nella pratica, infatti,  e importante sapere se un sistema  e o meno a banda limitata
e qual’ e la sua banda, per poter stabilire il campionamento pi u opportuno per le
immagini e, di conseguenza, per stabilire le caratteristiche del rivelatore.
Nel seguito si suppone che la PSF sia spazio-invariante e tale che
1. K(x) 0
2.
R
K(x)dx = 1 )
^
K(0) = 1
Le due condizioni implicano che la PSF  e integrabile e quindi, grazie al Teorema
di Riemann-Lebesgue, la funzione di trasferimento non solo  e limitata e continua
ma tende a zero per grandi frequenze cio e
j
^
K(!)j! 0 se j!j! +1:
Ne segue che un sistema di formazione di immagini  e un ltro di Fourier passa
basso (cio e le alte frequenze vengono soppresse). Inoltre la seconda condizione
implica che
Z
g(x)dx =
Z
f(x)dx
5
1 Perturbazioni delle immagini digitali
che  e detta condizione di conservazione dell’intensit a o di conservazione del usso.
Esempi di blurring sono:
 immagini mosse, che si formano quando un oggetto in movimento viene
ripreso da una camera ssa o quando un oggetto sso viene ripreso da una
camera in movimento; in entrambi i casi, durante il tempo di esposizione si
sovrappongono pi u immagini dello stesso ogetto, traslate nella direzione del
moto relativo;
 immagini sfocate, cio e non a fuoco;
 immagini con eetti di dirazione , che si formano soprattutto nelle acquisi-
zioni di microscopi e telescopi;
 immagini con blurring atmosferico, che sono dovute a turbolenze atmosferi-
che; il blurring atmosferico  e predominante nelle immagini prodotte dai gran-
di telescopi a terra ed  e stato recentemente superato tramite l’introduzione
dell’ottica adattiva.
1.2 Noise
Nel processo di formazione delle immagini digitali, la grandezza sica che rap-
presenta l’oggetto incognito f(x), che  e un segnale analogico, deve essere con-
vertito in un segnale discreto g(x), detto segnale digitale, che ne rappresenta la
distribuzione spaziale. Tale conversione viene solitamente eettuata mediante ar-
ray di sensori, ognuno dei quali fornisce una misura locale del valore di g(x).
Il segnale viene quindi prima trasformato da continuo a discreto (digitale), trami-
te la discretizzazione del prodotto (1.3). Si deve poi aggiungere un errore di tipo
casuale (quindi non noto): esso tiene conto del fatto che, per l’eetto di svariati
6
1.2 Noise
fattori, ogni sensore non fornisce un valore esatto di g, ma un suo valore appros-
simato.
Si vede prima la questione della discretizzazione e, per sempliticit a di notazione,
consideriamo solo il caso 2D; l’estensione al caso tridimensionale non presenta
dicolt a. Le immagini digitali sono discrete, quindi la funzione g(x) di variabili
continue  e sostituita da una matrice N N di valori dig(x) presi nei punti di un
reticolo regolare. Questi valori sono denotati con g
m;n
. Anche l’oggetto f
(0)
(x)
viene sostituito da una matriceN N di valori denotati conf
(0)
k;l
. A questo punto
la convoluzione integrale pu o essere sostituita da una convoluzione ciclica discre-
ta. Questa approssimazione  e ragionevole se la PSF  e signicativamente diversa
da zero su un dominio pi u piccolo di quello del dominio di partenza, cos  che gli
eetti della periodicit a siano trascurabili. Discretizzando la (1.3) nel caso delle
immagini 2D si ottiene
g
m;n
=
N1
X
k;l=0
K
mk;nl
f
k;l
(1.6)
o anche
g =K f: (1.7)
Occorre introdurre qui la notazione ‘vettorizzata’: si vede facilmente che il ri-
sultato della convoluzione ciclica (1.7) si pu o ottenere facilmente attraverso un
prodotto matrice per vettore. A tal scopo  e suciente trasformare le matrici in
vettori tramite un opportuno riordinamento degli elementi. Con l’ordine lessi-
cograco il vettore  e costruito accodando le righe l’una all’altra, per indice cre-
scente, e trasponendo il vettore riga nale cos  ottenuto; con l’ordinamento per
colonne (detto anche column-major o FORTRAN-like) vengono accodate le co-
lonne, per indice crescente. In questa tesi viene adottata la seconda convenzione
e si indica con ‘vect’ l’operatore corrispondente: dato un array M-dimensionale
V 2R
n
1
n
2
:::n
M
,v=vect(V )  e il vettore colonna con N = n
1
n
2
:::n
M
elementi.
7
1 Perturbazioni delle immagini digitali
In Matlab questa scelta  e vantaggiosa in quanto l’ordine degli elementi nel vetto-
re v corrisponde esattamente all’ordine dei corrispondenti valori nella memoria.
Inoltre, l’operatore vect non richiede alcuna implementazione particolare: il vet-
torev si ottiene infatti con l’operazione di indirizzazione degenere: [v=V(:).].
In questo contesto, sar a inoltre spesso utile riferirsi ad un determinato elemento
dell’array M-dimensionale V attraverso il suo indice lineare n invece che attra-
verso il vettore di indici (i
1
;i
2
;:::;i
M
) che identica la sua posizione in v. Con la
convezione qui adottata per l’operatore vect le relazioni fra le due indicizzazioni
sono le seguenti, dove si assume che tutti gli indici partano da 1:
n = (i
M
 1)n
1
n
2
 n
M1
+(i
M1
 1)n
1
n
2
 n
M2
+:::+(i
3
 1)n
1
n
2
+(i
2
 1)n
1
+i
1
(1.8)
i
e
=
 (n k
e
)=maxf1;m
e
g
 ; l =M;:::; 1 (1.9)
dove
m
e
=
 0 l = 1
P
l1
j=1
n
j
l> 1
e
k
e
=
 m
e+1
(i
e+1
 1) +k
e+1
1 l<M
1 l =M
In Matlab, la conversione dal multiindice (i
1
;:::;i
M
) all’indice linearen  e eettuata
dalla funzione sub2ind mentre la conversione opposta  e eettuata dalla funzione
ind2sub.
Si ha quindi che g(n) indica il valore esatto di g nel pixel n (od eventualmente il
suo integrale sul dominio del pixel) ed
f(n)
il valore esatto dell’oggetto nel pixel n.
Eettuando la discretizzazione si introduce gi a un’approssimazione rispetto all’og-
8
1.2 Noise
getto originale; la versione discreta della (1.3)  e data da
g(n) =
X
n
K(n n
0
)
 f(n
0
)
dove
g(n
1
;n
2
) =
N
X
n
1
=0
N
X
n
2
=0
K(n
1
 n
0
1
;n
2
 n
0
2
)
 f(n
0
1
;n
0
2
):
Inoltre si deve tener conto del fatto che il valore nel pixeln, fornito dal rivelatore,
non coincide esattamente cong(n) a causa dei fenomeni di carattere aleatorio che
vanno sotto il nome di noise. Pertanto, il valore misurato g(n)  e la realizzazione
di una variabile aleatoria  (n) il cui valore atteso  e pari a
E( (n)) = g(n)
.
Quantit a di noise
La quantit a di noise presente in un’immagine  e misurata dal rapporto segnale-
noise (signal-to-noise ratio, SNR), denito da
SNR = 10 log
10
 varianza dell’immagine
varianza del noise
 : (1.10)
Si misura in decibel (dB) e pu o essere stimato da
SNR = 20 log
10
 kgk
kwk
 (1.11)
con la normajjjj
1
oppurejjjj
2
. Se SNR  e dell’ordine di 40 o 50 dB, la com-
ponente di noise non  e visibile nell’immagine e gli eetti del blurring sono quelli
dominanti.
Sono presentati ora nel dettaglio i vari tipi di noise, che spesso compaiono simul-
taneamente.
9
1 Perturbazioni delle immagini digitali
Noise di lettura
Il noise di lettura (readout noise, RON )  e dovuto al sistema elettronico che
amplica il segnale in uscita ed  e di tipo addittivo; si pu o pertanto scrivere
g(n) = g(n) +w(n)
dove w(n)  e la realizzazione di una variabile aleatoria  che non dipende n e dal
segnaleg(n) n e dal pixel n. Ne segue quindi che la (1.3) e la sua versione discreta
(1.6) sono date da
g(x) =
Z
K(x x
0
)f(x
0
)dx
0
+w(x) (1.12)
e analogamente
g
m;n
=
N1
X
k;l=0
K
mk;nl
f
k;l
+w
m;n
: (1.13)
In genere il modello pi u usato per questo tipo di noise  e quello Gaussiano bianco,
secondo cui   e una variabile aleatoria Gaussiana con media 0 e varianza  :
P
 (w) =
1
p
2 e
 w
2
2 2
: (1.14)
Si sceglie questo modello perch e un graco della distorsione dei valori dei pixel
mostra una distribuzione normale del noise e perch e grazie al teorema del li-
mite centrale si ha che la somma di dierenti noise tende ad una distribuzione
gaussiana.
Noise fotonico
Se si trascura il noise di lettura, il valore misurato dell’immagine  e dato dal
numero di fotoni che, durante il tempo di esposizione, sono stati rilevati (noise
fotonico). g(n)  e modellizzato come la realizzazione di una variabile aleatoria
poissoniana con valore atteso g(n)
P
 (n)
(g(n)) =
e
 g(n)
g(n)
g(n)
g(n)!
(1.15)
10
1.2 Noise
dove g(n), proporzionale aT ,  e il valore atteso. Anche in questo caso si puo pen-
sare di scrivere g(n) nella forma additiva, ma bisogna ricordare che, a dierenza
del noise di lettura, w(n) non  e indipendente dal segnale.
Nel caso di conteggi alti, si pu o comunque approssimare con una gaussiana.
Noise sale e pepe
Il caso di noise sale e pepe  e qello in cui si hanno pixel sparsi molto luminosi
o oscuri, in cui l’immagine  e molto diversa per colore o intensit a rispetto ai pixel
circostanti. In genere questo tipo di noise coinvolge solo un piccolo numero di pi-
xel dell’immagine, che creano l’eetto di puntini bianchi nelle zone scure e puntini
neri nelle zone chiare.
Questo tipo di disturbo pu o essere causato da punti del rilevatore non pi u funzio-
nanti, da errori di conversioni da analogico a digitale o da errori nella trasmissione,
ecc.
11