Total Variation image denoising from Poisson data

Total Variation image denoising from Poisson
data
Split Bregman and Alternating Extragradient methods
Davide Taviani
22 luglio 2011

View Slide

Introduzione
• Ricostruzione di immagini degradate da rumore Poissoniano
• Nell’ambito del paradigma Bayesiano il problema si riduce alla
soluzione del problema di minimo vincolato
min
x≥0
φ(x) := φ0
(x) + βφ1
(x)
β parametro di regolarizzazione
φ0
(x) funzione di discrepanza dell’oggetto dai dati
φ1
(x) funzione di regolarizzazione

View Slide

Introduzione
• Regolarizzazione edge-preserving per preservare i contorni
delle aree presenti nell’immagine: scelta della Total Variation
come funzione di regolarizzazione (ben nota nel caso di rumore
Gaussiano)
φ1
(x) =
i
||(∇x)i
||2
• Metodi specializzati per il problema: metodi Split Bregman
• Formulazione primale-duale del problema ⇒ metodi di tipo
extragradiente

View Slide

Formazione dell’immagine ed errori
Immagine: segnale digitale ottenuto tramite uno strumento di
acquisizione.
Durante il processo di acquisizione si veriﬁcano errori:
• blurring: Fenomeno che pu`
o essere modellizzato in modo
deterministico
• rumore (noise): Errore di tipo aleatorio dovuto al processo di
memorizzazione dell’immagine

View Slide

Rumore
Due modelli sono comunemente usati per la descrizione del rumore:
• rumore Gaussiano (additivo): modello pi`
u diffuso per il
rumore, ampiamente usato per descrivere il rumore dovuto
all’agitazione termica degli elettroni in un conduttore (thermal
noise)
• rumore Poissoniano (rumore fotonico): modello per descrivere
il rumore presente in tutti quei processi in cui l’immagine si
ottiene tramite un conteggio di fotoni (microscopia ottica,
astronomia ottica/infrarossi, tomografia ad emissione, radiografia
digitale)

View Slide

Rumore Poissoniano
Probabilit`
a di ricevere k particelle in un dato intervallo di tempo T:
e−λλk
k!
k = 0, 1, 2, ...
λ valore atteso, proporzionale a T.
Rapporto segnale-rumore:
√
k.
⇓
Il rumore diventa importante nel caso di bassi conteggi (condizioni di
scarsa luminosit`
a).

View Slide

Massima verosimiglianza e KL
Seguendo l’approccio di Geman & Geman (1984), la stima di massima
verosimiglianza dell’immagine originale `
e data da ogni massimizzatore
della funzione di verosimiglianza:
Ly
(x) =
i
e−xi xi
fi
fi
!
Se ne considera il logaritmo negativo:
φ0
(x) =
i
fi
log
fi
xi
+ xi
− fi
Divergenza di Kullback-Leibler
Una soluzione banale a min
x
φ0
(x) `
e x = f (immagine rumorosa), non
soddisfacente.

View Slide

Regolarizzazione
Problema di ricostruzione di immagini: mal posto e mal
condizionato.
Per la regolarizzazione si usano ulteriori informazioni a priori sulla
natura dell’oggetto.
Approccio Bayesiano: si suppone che l’oggetto iniziale sia la
realizzazione di una variabile aleatoria X con distribuzione di
probabilit`
a
PX
(x) =
1
Z
e−βφ1(x)
• Z: costante di normalizzazione
• β: parametro di regolarizzazione
• φ1
(x): funzione di regolarizzazione

View Slide

Regolarizzazione
Possibili scelte per i parametri di regolarizzazione:
• β: un possibile approccio `
e il principio di discrepanza
recentemente proposto da Bertero et al. (2010)
• φ1
(x):
Penalizzazione Quadratica φ1
(x) =
||x||2
2
2
Laplaciano Quadratico φ1
(x) =
||∆x||2
2
2
Total Variation φ1
(x) =
i
||(∇x
)i
||2

View Slide

Formulazione del problema
Caso 2D: immagine con n pixel.
Se f `
e l’immagine rumorosa rilevata dal sistema di acquisizione, il
problema `
e formulato come:
min
x∈X
n
i=1
fi
log
fi
xi
+ xi
− fi
+ β
n
i=1
||Ai
x||2

View Slide

Formulazione del problema
Insieme ammissibile X:
X =



x ≥ η > 0 fi
> 0
x ≥ 0 fi
= 0
Ai
: matrice 2 × n con tutti gli elementi nulli tranne:
a1,i
= −1 a1,i+1
= 1
(Diﬀerenza tra un pixel e quello a destra)
a2,i
= −1 a1,i+t+1
= 1
(Diﬀerenza tra un pixel e quello sotto)

View Slide

Osservazioni sulla formulazione del problema
• φ0
(x) `
e convessa, coerciva e non lineare
• φ1
(x) `
e convessa e non ovunque differenziabile
⇓
• φ(x) `
e una funzione convessa, coerciva, non negativa e non
ovunque differenziabile
• la minimizzazione di φ(x) pu`
o comportare difficolt`
a sia di tipo
teorico che numerico (impossibilit`
a dell’uso di metodi per
funzioni differenziabili)
⇓
• sono stati appositamente progettati dei metodi per questo tipo di
problemi

View Slide

Iterazione di Bregman
Supponiamo di voler risolvere il seguente problema di ottimo
vincolato:
min
x
E(x) tale che Ax = b (1)
con E non diﬀerenziabile.
Usando una funzione di penalizzazione quadratica il problema si
riconduce a un problema di ottimo non vincolato:
min
x
E(x) +
1
2γ
||Ax − b||2
2
con γ costante positiva.

View Slide

L’idea `
e quella di sostituire, alla k-esima iterazione, alla parte non
diﬀerenziabile E la distanza di Bregman di E in x(k):
Dp
E
(x, x(k)) = E(x) − E(x(k) − p(k)
T
x − x(k)
dove p(k) denota il subgradiente di E in x(k).
L’iterazione di Bregman per questo problema `
e
x(k+1) = argmin
x
Dp
E
x, x(k) +
1
2γ
||Ax − b||2
2
= argmin
x
E(x) − p(k)T
x − x(k) +
1
2γ
||Ax − b||2
2
p(k+1) = p(k) −
1
γ
AT Ax(k+1) − b

View Slide

La garanzia che l’iterato x(k) giunga arbitrariamente vicino alla
soluzione del problema di minimo vincolato ci `
e dato dal seguente
teorema:
Teorema:
Sia E : Rn → R convessa; sia A : Rn → Rm lineare. Si consideri il
precedente algoritmo e si supponga che esista un iterato x∗ tale che
soddisfa Ax∗ = b.
Allora x∗ `
e una soluzione del problema di ottimo vincolato (1).

View Slide

p(k) ∈ sottospazio generato dalle colonne di AT
⇓
l’iterazione di Bregman si pu`
o riscrivere in maniera esplicita:
x(k+1) = argmin
x
E(x) +
1
2γ
Ax − b(k)
2
2
b(k+1) = b(k) + b − Ax(k+1)

View Slide

Metodi Split Bregman
Supponiamo di voler risolvere il seguente problema di minimo non
vincolato
min
x
{|Φ(x)| + H(x)}
con H e |Φ| funzioni convesse, e H e Φ differenziabili.
Una buona strategia `
e quella di isolare la parte differenziabile da
quella non differenziabile nel problema, introducendo il vincolo
artificiale Φ(x) = d.
In questo modo il problema diventa tale da poter applicare
l’iterazione di Bregman:
min
x,d
{|d| + H(x)} tale che d = Φ(x)

View Slide

L’iterazione di Bregman, in questo caso assume la seguente forma:
x(k+1), d(k+1) = argmin
x,d
|d| + H(x) +
1
2γ
d − Φ(x) − b(k)
2
2
(2)
b(k+1) = b(k) + Φ x(k+1) − d(k+1) (3)
L’espressione (2) pu`
o essere calcolata eﬃcientemente separando la
minimizzazione di x e d (da cui il termine split).

View Slide

Implementazione dell’iterazione Split Bregman:
while x(k) − x(k+1)
2
> tol
x(k+1) = argmin
x
H(x) + 1
2γ
d(k) − Φ(x) − b(k)
2
2
d(k+1) = argmin
d
|d| + 1
2γ
d − Φ x(k+1) − b(k)
2
2
b(k+1) = b(k) + Φ x(k+1) − d(k+1)
end
Per il calcolo di d(k+1) si pu`
o usare un cosiddetto operatore di
shrinkage.
La convergenza di questo metodo ci `
e assicurata dal fatto che `
e
possibile ricondurlo al metodo delle direzioni alternate dei
moltiplicatori, la cui convergenza `
e stata dimostrata (Gabay, 1983).

View Slide

Setzer et al. (2010) hanno proposto l’algoritmo PIDSplit+ per la
ricostruzione di immagini degradate da rumore Poissoniano, cio`
e un
implementazione della precedente iterazione di Bregman che consideri
le seguenti funzioni:
Divergenza di Kullback-Leibler (KL)
n
i=1
fi
log
fi
xi
+ xi
− fi
Total Variation (TV)
n
i=1
||Ai
x||2

View Slide

PIDSplit+
Si vuole risolvere il seguente problema di minimo:
min
x
n
i=1
fi
log
fi
xi
+ xi
− fi
+ β
n
i=1
||Ai
x||2
+ lx≥η
dove lx≥η
rappresenta la funzione indicatrice.
• Introduzione di tre vincoli, uno per la KL, uno per la Total
Variation e uno per la funzione indicatrice:
w1
= x w2
= Ax w3
= x
• Aggiunta al problema di minimo dei termini di penalizzazione, in
modo da ottenere:
min
x
n
i=1
fi log
fi
(w1)i
+ (w1)i − fi + β||w2||2 + lx≥η +
1
2γ
||x − w1||
2
2+
+ ||Ax − w2||
2
2 + ||x − w3||
2
2

View Slide

PIDSplit+
Ad ogni passo occore calcolare:
x(k+1) = 2I + AT A −1
w(k)
1
− b(k)
1
+ AT w(k)
2
− b(k)
2
+ w(k)
3
− b(k)
3
w(k+1)
1
=
1
2
b(k)
1
+ x(k+1) − γ + (b(k)
1
+ x(k+1) − γ)2 + 4γb
w(k+1)
2
= argmin
w2
β||w2
||2
+
1
2γ
b(k)
2
+ Ax(k+1) − w2
2
2
w(k+1)
3
=



b(k)
3
+ y(k+1) se b(k)
3
+ x(k+1) ≥ η
0 altrimenti
b(k+1)
1
= b(k)
1
+ x(k+1) − w(k+1)
1
b(k+1)
2
= b(k)
2
+ Ax(k+1) − w(k+1)
2
b(k+1)
3
= b(k)
3
+ x(k+1) − w(k+1)
3
(AT A) `
e una matrice circolante: la risoluzione di sistema in x(k+1) si pu`
o
calcolare eﬃcientemente con una Trasformata Discreta di Fourier.

View Slide

Formulazione primale-duale del problema
Si vuole risolvere il seguente problema di minimo:
min
x
i
fi
log
fi
xi
+ xi
− fi
+ β
i
||Ai
x||2
Dalla propriet`
a della norma in 2
, se a ∈ R2:
||a||2
= sup zT a con z ∈ R2 tale che ||z||2
≤ 1
otteniamo la riformulazione primale-duale del problema:
min
x
max
y
F(x, y)
con
F(x, y) =
i
fi
log
fi
xi
+ xi
− fi
+ βyT Ax

View Slide

Osservazioni sulla formulazione primale-duale
• Il problema di minimo `
e ora diventato un problema di punto
sella, dato che F `
e concava rispetto a y e convessa rispetto a x.
• La garanzia dell’esistenza di almeno un punto sella `
e data dal
teorema min-max (Rockafellar, 1970).
• Il problema pu`
o essere visto come una disequazione
variazionale (Variational Inequality Problem, VIP); infatti una
soluzione del problema primale-duale v∗ := (x∗, y∗) soddisfa
anche:
(v − v∗)T Φ(v∗) ≥ 0 ∀v ∈ X × Y
con
Φ(v) = Φ(x, y) =




∇
x
F(x, y)T
− ∇
y
F(x, y)T





View Slide

Metodi extragradiente
Per la soluzione di disequazioni variazionali una classe di metodi
efficaci sono i metodi extragradiente, caratterizzati da:
• proiezioni ad ogni passo
• necessit`
a di un parametro per la lunghezza del passo con:
un valore fisso che dipende da una stima a priori della costante di
Lipschitz del problema (potenzialmente difficile)
un valore aggiornato adattivamente ad ogni passo

View Slide

Alternating Extragradient Method
Ad ogni passo `
e necessario calcolare:
¯
y(k) = PY
y(k) + αk
∇y
F x(k), y(k)
x(k+1) = PX
x(k) + αk
∇x
F x(k), ¯
y(k)
y(k) = PY
y(k) + αk
∇y
F x(k+1), y(k)
La lunghezza del passo αk
, per assicurare convergenza del metodo,
deve rispettare: αk
∈ [ α , α ] con 0 < α < α e



1 − 2αk
Ak
− 2α2
k
B2
k
≥
1 − 2αk
Ck
≥

View Slide

Alternating Extragradient Method
con ∈ (0, 1) costante e, per il problema di ricostruzione di immagini
degradate da rumore Poissoniano:
Ak
=
f
x(k)
−
f
x(k+1)
2
||x(k+1) − x(k)||2
Bk
=
βAx(k+1) − βAx(k)
2
x(k+1) − x(k)
2
Ck
= 0

View Slide

Osservazioni sul metodo AEM
• Il metodo esegue aggiornamenti di tipo Gauss-Seidel sulla
variabile primale e duale.
• Ogni iterazione `
e costituita da tre proiezioni che, nel caso in
esame, sono particolarmente semplici per la struttura dei domini.
• Il parametro αk
pu`
o essere determinato in modo adattivo,
implementando una procedura di backtracking:
partenza da valore tentativo
veriﬁca della condizione di convergenza
in caso negativo si riduce il passo e si reitera la procedura

View Slide

Simulazioni numeriche
Confronto dell’eﬃcienza di AEM e PIDSplit+ tramite simulazioni
numeriche svolte su due tipi di problemi:
• Problemi simulati: immagini prive di rumore, esso viene
simulato tramite la funzione imnoise di Matlab
• Problemi reali: immagini degradate in partenza da rumore.
Per i problemi reali non si ha a disposizione una versione priva di
rumore: `
e stata usata la cosiddetta ground truth come immagine di
riferimento per valutare l’esito della ricostruzione ed `
e stata necessaria
l’implementazione di criteri di arresto.

View Slide

Problemi simulati: LCR
Figura: LCR-1 (256 × 256)

View Slide

Problemi simulati: Airplane
Figura: Airplane (256 × 256)

View Slide

Problemi simulati: DR
Figura: DR (512 × 512)

View Slide

Problemi simulati
• Per β sono stati usati valori noti per fornire buoni risultati (R.
Zanella et al., 2009)
• Per il parametro γ nell’algoritmo PIDSplit+ sono stati usati i
seguenti valori: 50
β
, 5
β
, 1
β
, 0.5
β
• Entrambi gli algoritmi sono stati eseguiti con 3000 iterazioni
• `
E stata calcolata una soluzione “ideale” usando AEM con 10000
iterazioni

View Slide

Introduzione di due tipi di errore:
(iter) =
x(iter) − ˜
x
2
||˜
x||
2
errore rispetto alla soluzione ideale ˜
x
e(iter) =
x(iter) − x∗
2
||x∗||
2
errore di ricostruzione

View Slide

Prima serie di test:
• ﬁssate alcune soglie per ( 10−4 , 10−5 , 5 · 10−6 )
• misurate le iterazioni il tempo necessari a raggiungere , insieme
all’errore di ricostruzione ottenuto
• misurato il tempo necessario a compiere 3000 iterazioni e l’errore
di ricostruzione ottenuto

View Slide

Andamento di rispetto al tempo
0 20 40 60 80 100 120 140 160
10−5
10−4
10−3
10−2
AEM
PID50
PID5
PID1
PID05
0 20 40 60 80 100 120 140 160
10−7
10−6
10−5
10−4
10−3
10−2
AEM
PID50
PID5
PID1
PID05
LCR-1 Airplane
0 100 200 300 400 500 600
10−5
10−4
10−3
10−2
10−1 AEM
PID50
PID5
PID1
PID05
DR

View Slide

Seconda serie di test:
• misurato il numero di iterazioni compiute dopo 5 secondi insieme
all’errore di ricostruzione
• comparate le immagini ottenute dopo 5 secondi

View Slide

Andamento di e nei primi 5 secondi
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
0.03
0.04
0.05
0.06
0.07
0.08
0.09
AEM
PID50
PID5
PID1
PID05
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
0.02
0.025
0.03
0.035
0.04
0.045
0.05
0.055
0.06
0.065
0.07
AEM
PID50
PID5
PID1
PID05
LCR-1 Airplane
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
0.04
0.06
0.08
0.1
0.12
0.14
0.16
AEM
PID50
PID5
PID1
PID05
DR

View Slide

Riga 128 di LCR-10 dopo 5 secondi
0 50 100 150 200 250
200
400
600
800
1000
1200
1400
1600
1800
2000
0 50 100 150 200 250
200
400
600
800
1000
1200
1400
1600
1800
2000
AEM PID5
PIDSplit+ con γ = 5
β

View Slide

Conclusioni
• Forte inﬂuenza di γ nella velocit`
a di convergenza di PIDSplit+: i
valori con cui si ha una velocit`
a iniziale elevata non sono
soddisfacenti asintoticamente, e viceversa.
• Posizionamento intermedio di AEM rispetto al comportamento
migliore e peggiore di PIDSplit+, sia per la velocit`
a iniziale che
quella asintotica.
• Minor costo computazionale di ogni iterazione di AEM rispetto a
ogni iterazione PIDSplit+.
• Assenza di parametri forniti dall’utente per AEM.
• La procedura di backtracking per αk
viene invocata un numero
limitato di volte, soprattutto nelle iterazioni iniziali.
• Una scelte oculata del parametro γ in PIDSplit+ pu`
o fornire
buoni risultati anche pi`
u velocemente di AEM.

View Slide

Total Variation image denoising from Poisson data

Total Variation image denoising from Poisson data

Davide Taviani

More Decks by Davide Taviani

Other Decks in Science

Featured

Transcript