Total Variation image denoising from Poisson data

8cad6d8aa26abc031f3c5c38a80fd06e?s=47 Davide Taviani
January 04, 2012

Total Variation image denoising from Poisson data

8cad6d8aa26abc031f3c5c38a80fd06e?s=128

Davide Taviani

January 04, 2012
Tweet

Transcript

  1. 1.

    Total Variation image denoising from Poisson data Split Bregman and

    Alternating Extragradient methods Davide Taviani 22 luglio 2011
  2. 2.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    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 ≥
  27. 27.

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

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

    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.
  30. 33.

    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
  31. 34.

    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
  32. 35.

    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
  33. 36.

    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
  34. 37.

    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
  35. 38.

    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
  36. 39.

    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 β
  37. 40.

    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.