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
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
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: • 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)
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).
(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.
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
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
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
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)
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
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.
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
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).
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)
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).
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).
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
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
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
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
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
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
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
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.
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
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
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.
Davide Taviani
helios
iamctodd
reverentgeek
rmw
skipperchong
carmenhchung
keathley
productmarketing
clairelew
thoeni
