$30 off During Our Annual Pro Sale. View Details »

Total Variation image denoising from Poisson data

Davide Taviani
January 04, 2012

Total Variation image denoising from Poisson data

Davide Taviani

January 04, 2012
Tweet

More Decks by Davide Taviani

Other Decks in Science

Transcript

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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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)

    View Slide

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

    View Slide

  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.

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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)

    View Slide

  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

    View Slide

  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

    ||Ax − b||2
    2
    con γ costante positiva.

    View Slide

  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

    ||Ax − b||2
    2
    = argmin
    x
    E(x) − p(k)T
    x − x(k) +
    1

    ||Ax − b||2
    2
    p(k+1) = p(k) −
    1
    γ
    AT Ax(k+1) − b

    View Slide

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

    View Slide

  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

    Ax − b(k)
    2
    2
    b(k+1) = b(k) + b − Ax(k+1)

    View Slide

  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)

    View Slide

  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

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

    View Slide

  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

    d(k) − Φ(x) − b(k)
    2
    2
    d(k+1) = argmin
    d
    |d| + 1

    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

  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

    View Slide

  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

    ||x − w1||
    2
    2+
    + ||Ax − w2||
    2
    2 + ||x − w3||
    2
    2

    View Slide

  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

    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.

    View Slide

  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

    View Slide

  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




    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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.

    View Slide

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

    View Slide

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

    View Slide

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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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

    View Slide

  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
    β

    View Slide

  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.

    View Slide