Problemi Inversi Malposti - dm.unibo.itmontelau/html/Problemi Inversi Malposti.pdf · Problemi...

Serena Papi Esercitazione laboratorio -Matematica computazionale 2008/2009 1

Problemi Inversi Malposti

Consideriamo il sistema lineare

A x = b

(1)

dove la matrice A è nxn, ma è molto mal condizionata. Il problema (1) viene

detto problema inverso mal posto. Pertanto la soluzione del sistema (1) è

estremamente sensibile a perturbazioni, anche relativamente piccole, sul termine

noto b.

Consideriamone la versione con termine noto perturbato εε += bb ; avremo

allora che il sistema da risolvere è

εε bAx = .

Approccio SVD alla Soluzione Poiché siamo in presenza anche di rumore occorre filtrare la soluzione. Sia A una matrice nxn si definisce SingularValueDecomposition di A

A = UΣVT dove UTU = VTV = In Σ = diag(σ1, …, σn) con

σ1 ≥ σ2 … σn ≥ 0 Tutti i valori singolari decadono gradualmente e l’indice di condizionamento di A

cond(A)= σ1/ σn


è molto elevato perché la matrice A è molto mal condizionata. L’immediata soluzione che ne deriva è

xε = A-1 bε = i

N

i i

T

i vbu

∑=1 σ

ε =

=A-1 b + A-1 ε = i

N

i i

T

i vbu

∑=1 σ

+ i

N

i i

T

i vu

∑=1 σ

ε = x + i

N

i i

T

i vu

∑=1 σ

ε

ed è completamente dominata dall’inversione della componente di rumore A-1ε. Vediamo un esempio su un segnale monodimensionale per capire la mal

posizione del problema.

Si consideri il seguente segnale lungo n=512: >> x=load( 'segnale_box.dat' ); >> figure >> plot(x);

0 100 200 300 400 500 6000

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Si consideri ora una matrice A 512x512 che applicata al segnale ne produce una versione distorta


>> load -mat 'matrice_conv_mono.dat' >> b=A*x'; >> figure >> plot(b);

0 100 200 300 400 500 6000

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

>> be=addnoiseG(b,0,0.00001); >> figure >> plot(be);

0 100 200 300 400 500 600-0.2

0

0.2

0.4

0.6

0.8

1

1.2

La matrice A di distorsione è una matrice di Toeplitz molto sparsa generata da una stringa di valori che rappresentano il campionamento di una gaussiana sul suo supporto


Vediamo le soluzioni ottenute risolvendo il sistema applicando 1. Soluzione mediante QR: >>[Q,R] = qr(A); >> be_new = Q’*be; >> x_sol= R\be_new;

0 100 200 300 400 500 600-80

-60

-40

-20

0

20

40

60

80


2. Approccio mediante SVD [U,S,V]=svd(A); x_sol=zeros(N,1); for i=1:N x_sol=x_sol+((U(:,i)'*be)/S(i,i))*V(:,i); end

0 100 200 300 400 500 600-60

-40

-20

0

20

40

60

Vediamo come decadono I valori singolari di A >> figure >> plot(diag(S))

0 100 200 300 400 500 6000

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1


Metodi di Regolarizzazione 1. Approccio TSVD alla Soluzione TSVD: Un semplice approccio per ridurre il rumore in ricostruzione è quello di

scartare tutte le componenti della SVD che sono dominate dal rumore.

Questo porta al soluzione mediante TSVD:

xk = .,1

Nkvbu

i

k

i i

Ti <∑

= σ

Nonostante la sua semplicità questo metodo può lavorare abbastanza bene se si

riesce a determinare il giusto parametro k di troncamento.

Tornando all’esempio la soluzione che si ottiene è

x_sol=zeros(N,1); %inserire k il numero dei valori singolari da considerare k=128; for i=1:k x_sol=x_sol+((U(:,i)'*ye)/S(i,i))*V(:,i); end figure plot(x_sol)


k=100

0 100 200 300 400 500 600-0.2

0

0.2

0.4

0.6

0.8

1

1.2

k=128

0 100 200 300 400 500 600-0.2

0

0.2

0.4

0.6

0.8

1

1.2


k=180

0 100 200 300 400 500 600-0.2

0

0.2

0.4

0.6

0.8

1

1.2

K=250

0 100 200 300 400 500 600-0.2

0

0.2

0.4

0.6

0.8

1

1.2


2.Metodi di filtraggio spettrale La TSVD è un esempio di una classe generale di metodi, noti come metodi di

filtraggio spettrale, in cui la soluzione è

xfilt = i

N

i i

Ti

i vbu

∑=1 σ

φ ,

con filtri φi;

Tali filtri sono scelti in modo tale che φi∼1 per i grandi valori singolari e φi∼0

per i valori singolari piccoli.

Infatti nella TSVD si ha:

+====

Nki

ki

i

i

,...,10

,...,11

φφ

Dove il parametro k è noto come parametro di troncamento e determina il

numero di componenti della SVD da mantenere nella soluzione regolarizzata.

Regolarizzazione di Tikhonov

In questo metodo si definisce il fattore di filtro come

Nii

i

i ,...,122

2

=+

=λσ

σφ

con λ>0 parametro di regolarizzazione.


Quando

• σi >>λ si ha φi∼1

• σi =<λ si ha φi∼σi2/λ2

Il parametro λ nel metodo di Tikhonov e il parametro k nella TSVD agiscono

allo stesso modo ovvero controllano le componenti della SVD da filtrare.

Determinazione del parametro di filtraggio

Prendiamo in considerazione 3 importati criteri di determinazione del parametro

• Il principio di discrepanza

• Il metodo della cross-validation

• Il criterio della L-curva.

Si consideri la norma della soluzione filtrata

2

1

2

2 ∑=

=

N

i i

Ti

ifilt

bux

σφ ,

e la norma del residuo

( )2

1

2

2)1(∑

=

−=−N

i

Tiifilt buAxb φ .

Nel metodo TSVD si ha

• La norma della soluzione è una funzione di k monotonicamente nondecrescente

• La norma del residuo è invece monotonicamente

noncrescente. Nel metodo di Tikhonov al contrario si ha

• La norma della soluzione è una funzione di λ monotonicamente noncrescente


• La norma del residuo è invece monotonicamente nondecrescente.

Principio di Discrepanza

Tale principio richiede una buona stima della deviazione standard δ del rumore.

L’idea che sta alla base di tale principio è quella di scegliere il parametro di

regolarizzazione in modo che la norma del residuo sia approssimativamente pari

alla deviazione standard del rumore:

δτ=−2

2filtAxb

con τ>1 valore reale.

Data la SVD e i filtri UTb il costo di G(k) è di 2N moltiplicazioni e 2N addizioni

per ogni calcolo del residuo.

Il principio di discrepanza è convergente anche quando il rumore tende a zero

pero’ richiede informazioni che spesso non sono disponibili.

Tornando al nostro esempio il risultato che si ottiene applicando la

regolarizazione di Tikhonov in cui il parametro di regolarizzazione è

determinato mediante il principio della discrepanza è mostrato di seguito:


0 100 200 300 400 500 600-0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Il codice per la risoluzione è “approssimazione.m”


Applicazioni:

Image Restoration

X B immagine originale immagine blurred

(distorta) n x n n x n

Il modello di distorsione è un modello lineare

Ax = b Per poterlo rappresentare occorre procedere ad una vettorizzazione delle

immagini.

La vettorizzazione per colonne delle matrici-immagini X e B nei vettori x e b (di

dimensione n2 x 1) può essere facilmente realizzata in Matlab:

x = X(:);

X = reshape(x,n,n);

PSF


Modello Lineare: Esiste una matrice A di dimensioni N x N (N=n2) tale che

b = Ax, dove

A è la matrice discretizzazione della PSF (spesso identificata come matrice di blurring o distorsione)

B e X sono le matrici relative alle immagini b e x le loro vettorizzazioni.

Perciò dobbiamo risolvere un problema inverso, ovvero determinare x, noti b e A

Ax = b

Il problema dell’image restoration (deblurring) consiste nel risolvere questo sistema. La difficoltà è dovuta a differenti motivi:

• la matrice A è di grandi dimensioni • la matrice A è molto mal condizionata

• la matrice A potrebbe essere un modello di blurring non esatto


La funzione di distorsione (Point Spread

Function) Le immagini distorte nascono dai processi meccanici e fisici di acquisizione.

Un esempio tipico sono le immagini sfuocate prodotte da limitazioni del

sistema ottico oppure immagini mosse dovute ad un movimento dell’oggetto

ripreso.

La matrice A rappresenta il blurring che applicato all’immagine originale x

ne ha prodotto la versione distorta b.

Si vuole dare una caratterizzazione della matrice A di blurring.

Sorgente Puntuale Sorgente di dispersione puntuale

La PSF solitamente è locale: questo significa che la zona con elementi non

nulli è piccola.

E’ possibile quindi memorizzare e salvare in maniera sparsa.


Modelli di PSF Horizontal motion blur Vertical motion blur

Out of focus Atmospheric turbulence blur


• Horizontal motion blur e Vertical motion blur

• Out of focus

• Atmospheric turbulence blur


IN MATLAB: I = imread('peppers.png');

I = I(60+[1:256],222+[1:256],:); % crop the image

figure; imshow(I);

LEN = 31; THETA = 11;

PSF = fspecial('motion',LEN,THETA); % create PSF

Blurred = imfilter(I,PSF,'circular','conv');

figure; imshow(Blurred);


Il problema dei bordi La linea gialla delimita i bordi dell’immagine blurred B.

La PSF va oltre i bordi dell’immagine, e questo causa nell’immagine B una

perdita dell’informazione dell’immagine originale.

Osservando la figura si puo’ notare che, se si considera un pixel bij di bordo

nell’immagine blurred questo è ottenuto come somma pesata del pixel xij e di

quelli del suo intorno, alcuni dei quali possono esere fuori dal campo di vista.

Un buon modello di image deblurring deve tener conto del problema ai bordi

altrimenti la ricostruzione potrà contenere alcuni artefatti spiacevoli sui bordi.

La tecnica più usata per trattare con la mancanza di informazione è quella di

imporre alcune condizioni ai bordi.

Questo significa che definiamo delle condizioni per determinare Xext a partire da

X.

1. zero boundary conditions

Si assume che l’immagine esatta sia nera (zero) fuori dai bordi.

Xext =

000

0X0

000

0 = sottomatrice nulla


Tale scelta è buona quando l’immagine esatta è per lo più nulla ai bordi come spesso accade per le immagini astronomiche. Sfortunatamente questa scelta ha un brutto effetto su quelle immagini che invece non sono nulle vicino ai bordi. Infatti spesso appaiono alcuni artefatti neri ai bordi oppure degli anelli ai bordi causati dalla grossa differenza tra i valori dentro e fuori dal bordo.

Le due condizioni ai bordi che seguono impongono un modello più realistico sul comportamento delle immagini ai bordi facendo uso solo delle informazioni disponibili, ovvero delle informazioni dell’ immagine vicino al bordo.


2. periodic boundary conditions

L’immagine si ripete, come se fosse un segnale periodico, in

tutte le direzioni.

Xext =

XXX

X XX

XXX

Si ottiene il seguente effetto:


3. reflexive boundary condition

L’immagine fuori dalla scena è l’immagine specchiata dell’immagine interna.

Xext =

xudx

lrlr

xudx

XXX

X XX

XXX

dove, utilizzando la sintassi del matlab, Xx=fliplr(X); X ud=flipud(X); Xx=fliplr(X ud);

Serena Papi Esercitazione laboratorio -Matematica computazionale 2008/2009

23

Il rumore Oltre alla distorsione le immagini osservate solo solitamente presentano anche il rumore. Il rumore dipende da differenti sorgenti e può essere lineare, nonlineare, moltiplicativo o additivo. Solitamente si considera un modello di rumore additivo di tipo gaussiano che caratterizza i sistemi di conversione analogico-digitale. La funzione Matlab che realizza il rumore gaussiano è ‘randn’ E = 0.01*randn(n,n); E è una matrice n x n che contiene gli elementi di rumore gaussiano con media 0 e standard deviation 0.01. Perciò il modello finale risulta:

Bnoisy = Bext + E

Il corrispondente modello di deblurring ora diviene B = (A0 + ABC)X + E Nell’image processing toolbox vi è la funzione imnoise che data un’immagine gli applica rumore permettendo di scegliere il tipo di rumore e la quantità: Bnoisy =imnoise(B, type, parameters); Bnoisy =imnoise(B, ‘gaussian’, mean, variance);


24

1.1 TSVD

Codice Matlab: b=double(B(:)); % B matrice dell’immagine blurred e rumore [U,E,V]=svds(A,1000); Ut=U'; Y=Ut*b; C=(Y')/E; x=V*(C'); X=reshape(x,n,n); X=centra(X,estensione);

Immagine originale Immagine blurred Immag. blurred + rumore

Immagine filtrata con k=100 Immagine filtrata con k=1000


25

Il bordo nero è dovuto alla zero boundary condition. Sotto viene mostrato lo stesso esempio ma considerando la reflexive boundary condition:

Problemi Inversi Malposti - dm.unibo.itmontelau/html/Problemi Inversi Malposti.pdf · Problemi...

Documents

Transcript of Problemi Inversi Malposti - dm.unibo.itmontelau/html/Problemi Inversi Malposti.pdf · Problemi...