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

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 

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 

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

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) 

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). 

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. 

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 

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 

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 

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 (Differenza tra un pixel e quello a destra) a2,i = −1 a1,i+t+1 = 1 (Differenza tra un pixel e quello sotto) 

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 

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 differenziabile. 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. 

Iterazione di Bregman L’idea ` e quella di sostituire, alla k-esima iterazione, alla parte non differenziabile 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 

Iterazione di Bregman 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). 

Iterazione di Bregman 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) 

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) 

Metodi Split Bregman 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 efficientemente separando la minimizzazione di x e d (da cui il termine split). 

Metodi Split Bregman 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). 

Metodi Split Bregman 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 

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 

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 efficientemente con una Trasformata Discreta di Fourier. 

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 

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     

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 

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 ≥ 

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 

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 verifica della condizione di convergenza in caso negativo si riduce il passo e si reitera la procedura 

Simulazioni numeriche Confronto dell’efficienza 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. 

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

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

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

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 

Simulazioni numeriche 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 

Simulazioni numeriche Prima serie di test: • fissate 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 

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 

Simulazioni numeriche 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 

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 

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 β 

Conclusioni • Forte influenza 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.