La matematica di diffuse or sharpen

Questo articolo documenta la matematica attualmente implementata nel modulo diffuse or sharpen di Ansel, così come si trova in src/iop/diffuse.c, src/common/bspline.h, data/kernels/diffuse.cl e data/kernels/bspline.cl. Non è una guida all’uso. È una ricostruzione del modello numerico a partire dal codice sorgente, con le affermazioni scientifiche ricondotte ai riferimenti citati nei commenti del codice.¹

Sommario

Il modulo opera su una decomposizione multiscala à-trous a piena risoluzione dell’immagine, costruita da sfocature con B-spline cardinali. A ogni scala, Ansel calcola quattro operatori anisotropi del secondo ordine separatamente sulle bande a bassa frequenza e ad alta frequenza, ne regolarizza l’azione tramite un’energia locale della banda ad alta frequenza normalizzata rispetto alla scala, restringe opzionalmente l’aggiornamento a una maschera binaria di inpainting e ricostruisce l’immagine dal grossolano al fine.² Gli operatori spaziali discreti sono stencil a differenze centrate 3x3 campionati sul sotto-reticolo sparso à-trous della scala corrente, per cui il supporto fisico più ampio deriva dal passo $2^s$ della scala wavelet piuttosto che da un kernel PDE più grande.³¹

Motivazione

L’obiettivo del modulo diffuse era simulare la fuoriuscita di pigmento tipo acquerello da un’immagine verso i suoi bordi, come questa :

Si è scoperto che il semplice cambio di segno dell’aggiornamento dell’equazione alle derivate parziali di diffusione poteva in realtà aumentare la nitidezza dell’immagine. Avevo lavorato per anni sulla deconvoluzione cieca quando iniziai a lavorare su un modello di diffusione, e non riuscii mai a far funzionare nemmeno il metodo allo stato dell’arte al di fuori dei bei casi da manuale (basso ISO, bassa saturazione, sfocatura uniforme). Ma questo naturalmente richiedeva un parametro di regolarizzazione per evitare la divergenza.

È così che nacque il modulo diffuse or sharpen, come framework generico per la simulazione di ogni sorta di fenomeno diffusivo o contro-diffusivo.

Modello continuo

In sostanza, il modulo implementa una famiglia di equazioni di diffusione anisotropa della forma

$$ \partial_t u = \nabla \cdot (\mathbf{A} \nabla u), $$

dove $u$ è l’immagine e $\mathbf{A}$ è un tensore di diffusione simmetrico e definito positivo che orienta la levigatura sia lungo le isofote, sia lungo il gradiente, sia in modo isotropo.²⁴

L’ispirazione originale nel codice è il modello di inpainting di Qin et al., che accoppiano il ripristino della struttura e della texture dell’immagine tramite il trasferimento anisotropo del calore.² Ansel mantiene quello spirito, ma applica la PDE in un dominio wavelet, separatamente sulle bande a bassa frequenza e ad alta frequenza, ed espone quattro coefficienti di trasporto controllati dall’utente.

Gradienti discreti e anisotropia

Per ogni pixel e ogni canale, il modulo estrae un intorno 3x3 e valuta le derivate prime centrate

$$ \begin{align} g_x &= \frac{u(i+1,j) - u(i-1,j)}{2} \\ g_y &= \frac{u(i,j+1) - u(i,j-1)}{2} \end{align} $$

Questo è il gradiente a differenze centrate standard usato nei metodi discreti di scale-space e di diffusione.³

L’orientamento locale del gradiente è quindi :

$$ \begin{align} \cos\theta = \frac{g_x}{\sqrt{g_x^2 + g_y^2}} \\ \sin\theta = \frac{g_y}{\sqrt{g_x^2 + g_y^2}} \end{align} $$

con il consueto ripiego su $(1,0)$ quando la magnitudine del gradiente si annulla.¹ Nota che eviteremo di calcolare il costoso angolo $\theta$ come la funzione $\arctan2$ delle componenti $(x, y)$, poiché non verrà mai usato direttamente nel seguito.

L’intensità dell’anisotropia è convertita dal parametro utente $a$ in un coefficiente positivo

$$ \alpha = a^2, $$

mentre il segno di $a$ seleziona la modalità:

• $a = 0$: diffusione isotropa,
• $a > 0$: diffusione allineata alla direzione dell’isofota,
• $a < 0$: diffusione allineata alla direzione del gradiente.¹

Il termine di smorzamento è quindi

$$ c^2 = \exp(-\alpha \lVert \nabla u \rVert), $$

che corrisponde alla diffusione anisotropa di Qin et al. per l’inpainting.²

Per la diffusione allineata alle isofote, il tensore scritto nella base locale dell’immagine è :

$$ \mathbf{A}_{\perp} = \begin{bmatrix} \cos^2\theta + c^2 \sin^2\theta & (c^2 - 1)\cos\theta\sin\theta \\ (c^2 - 1)\cos\theta\sin\theta & c^2 \cos^2\theta + \sin^2\theta \end{bmatrix}. $$

Per la diffusione allineata al gradiente, Ansel usa la forma invertita :

$$ \mathbf{A}_{\parallel} = \begin{bmatrix} c^2 \cos^2\theta + \sin^2\theta & (1 - c^2)\cos\theta\sin\theta \\ (1 - c^2)\cos\theta\sin\theta & \cos^2\theta + c^2 \sin^2\theta \end{bmatrix}. $$

Gli stencil di diffusione 3x3

A seconda dei parametri utente, valuteremo per ogni pixel lo stencil $\mathbf{K(A)}$ usando uno dei tensori $\mathbf{A}$ sopra.

Dato il tensore simmetrico $\mathbf{A}$ :

$$ \mathbf{A} = \begin{bmatrix} a_{11} & a_{12} \\ a_{12} & a_{22} \end{bmatrix}, $$

costruiamo il kernel Laplaciano anisotropo, ruotato e discreto $\nabla \cdot (\mathbf{A}\nabla u)$ come lo stencil 3x3 :

$$ \mathbf{K}(\mathbf{A}) = \begin{bmatrix} \frac{a_{12}}{2} & a_{22} & -\frac{a_{12}}{2} \\ a_{11} & -2(a_{11}+a_{22}) & a_{11} \\ -\frac{a_{12}}{2} & a_{22} & \frac{a_{12}}{2} \end{bmatrix}, $$

a meno della convenzione di segno menzionata nel codice sorgente per i termini fuori diagonale.²⁴¹

Nel caso isotropo si riduce a una costante :

$$ \mathbf{K}_{\text{iso}} = \begin{bmatrix} \tfrac14 & \tfrac12 & \tfrac14 \\ \tfrac12 & -3 & \tfrac12 \\ \tfrac14 & \tfrac12 & \tfrac14 \end{bmatrix}, $$

che è il classico Laplaciano isotropo in stile Oono-Puri. Il suo interesse sta nel comportamento rotazionale: rispetto al Laplaciano a 5 punti, l’errore angolare è ridotto, il che è rilevante quando la diffusione non deve privilegiare gli assi dell’immagine.⁵⁶ Questo errore angolare è stato valutato rispetto a kernel simili in Rotation-invariant Laplacian for 2D grids .

Per comprendere l’effetto della direzione nella diffusione, partiamo da un disco rumoroso e diffondiamolo con Ansel :

Un paio di cose vale la pena menzionare qui :

• Solo la sfocatura isotropa produce un disco acromatico, le altre trattengono un po’ di rumore cromatico su larga scala, perché il rumore cromatico crea variazioni locali del gradiente.
• La diffusione anisotropa è perfetta vicino alle posizioni nord/est/sud/ovest del cerchio, cioè quando l’angolo del gradiente è perfettamente allineato alla griglia dei pixel (0° o 90°). Nelle posizioni intermedie, vediamo discrepanze riguardo a come viene trattata la direzione, dovute ai limiti numerici di uno stencil quadrato 3×3 ruotato.

La piramide B-spline à-trous

L’analisi multiscala è costruita a partire dal filtro B-spline cardinale separabile a 5 tap

$$ h_0 = \frac{1}{16}[1,4,6,4,1]. $$

Questo filtro è un’approssimazione gaussiana compatta, motivo per cui compare naturalmente nei metodi scale-space basati su spline.⁷ Riutilizziamo qui il framework introdotto da Johannes Hanika per il modulo contrast equalizer.⁸

L’implementazione attuale usa la storica scala à-trous non decimata:

1. Il primo passa-basso a piena risoluzione è

   $$ G_0 = h_0 * u, $$

2. La banda di dettaglio più fine è

   $$ H_0 = u - G_0. $$

3. I livelli più grossolani mantengono la stessa risoluzione dell’immagine e si limitano ad allargare il passo della sfocatura di $2^s$:

   $$ G_s = h_{2^s} * G_{s-1}, \qquad s > 0. $$

   ovvero a $s=1$ i tap $h_{2}$ sono valutati ogni pixel successivo, a $s=2$ i tap $h_4$ sono valutati ogni 4 pixel, e così via…

4. La banda memorizzata alla scala $s$ è la differenza tra due passa-basso successivi:

   $$ H_s = G_{s-1} - G_s, \qquad s > 0. $$

Così ogni banda vive sulla griglia originale dell’immagine (nessuna decimazione), e risolviamo la PDE dal fine al grossolano in parallelo, accumulando la soluzione scala per scala nel buffer di output finché non aggiungiamo il residuo finale.¹ Questo schema evita di impilare la correzione delle scale fini sopra le scale grossolane, come è consuetudine con le piramidi gaussiane o i solutori multigrid, e ha mostrato empiricamente una migliore stabilità nel problema inverso della ricostruzione della nitidezza. Ciò equivale a gestire l’energia di ciascuna banda separatamente, in un modo che ricorda gli equalizzatori audio HiFi, ma all’interno di un framework spaziale 2D.

Scala del kernel, scala della banda e inviluppo di scala della GUI

Per studiare le proprietà dello schema wavelet, useremo la proprietà della B-spline cardinale di essere un’approssimazione di filtro gaussiano. Il parametro gaussiano $\sigma$ controlla come si compongono i kernel di sfocatura:

$$ G(\sigma_1) * G(\sigma_2) = G\left({\sqrt{\sigma_1^2 + \sigma_2^2}}\right) $$

per cui le varianze $\sigma^2$ dei kernel gaussiani si sommano sotto convoluzione.³⁷ Al contrario, la varianza di un segnale filtrato $X$ :

$$ \operatorname{Var}[g_\sigma * X], $$

dipende dallo spettro di $X$ e in generale non è uguale a $\sigma^2$. Per disambiguare entrambe le varianze, nel resto di questo articolo $\sigma$ si riferisce solo al parametro del kernel gaussiano, mai alla radice quadrata della varianza del segnale.

Per evitare ambiguità, aiuta anche separare i livelli passa-basso dalle bande dei livelli:

• $G_0 = h_1 * u$ è il primo passa-basso a piena risoluzione sull’immagine di input $u$,
• $G_s = h_{2^s} * G_{s-1}$ per $s \ge 1$,
• $H_0 = u - G_0$,
• $H_s = G_{s-1} - G_s$ per $s \ge 1$.

Quindi:

• l’indice passa-basso $s$ si riferisce al livello di sfocatura $G_s$,
• l’indice di banda $s$ si riferisce alla banda di dettaglio $H_s$ situata tra $G_{s-1}$ e $G_s$,
• ogni livello rimane campionato sulla griglia originale dell’immagine.¹

Il filtro B-spline approssima al meglio il kernel gaussiano equivalente di parametro :

$$ \sigma_B \approx 1.055365. $$

Poiché le varianze dei kernel gaussiani si sommano sotto convoluzione, il raggio di sfocatura effettivo di $G_s$ segue direttamente dalla scala di analisi à-trous. Pertanto la sequenza cumulativa dei passa-basso è

$$ G_0; G_1; G_2; G_3; \dots \quad \Longleftrightarrow \quad \sigma_{G_0}; \sigma_{G_1}; \sigma_{G_2}; \sigma_{G_3}; \dots $$

con

$$ \sigma_{G,0} = \sigma_B,\qquad \sigma_{G,1} = \sqrt{5}\sigma_B,\qquad \sigma_{G,2} = \sqrt{21}\sigma_B,\dots $$

e, in generale,

$$ \sigma_{G,s}^2 = \sum_{k=0}^{s} \sigma_B^2 4^k = \sigma_B^2 \frac{4^{s+1} - 1}{3}. $$

Equivalentemente,

$$ \sigma_{G,s} = \sigma_B \sqrt{\frac{4^{s+1} - 1}{3}}, $$

Il codice usa $\sigma_{G_s}$ per decidere quante scale sono necessarie per soddisfare il parametro radius richiesto dall’utente e per pesare ogni banda attorno al raggio centrale selezionato dall’utente:

$$ w_s = \exp\left( - \frac{(z\sigma_{G,s} - r_c)^2}{r_w^2} \right), $$

dove $z$ è il livello di zoom della camera oscura (rispetto al raw a piena risoluzione), $r_c$ è radius_center nella GUI, e $r_w$ è radius.¹ Impostando questi valori, gli utenti definiscono il decadimento di una sorta di filtro di “larghezza di banda” discretizzato, centrato su una frequenza arbitraria. In questo senso, $r_w$ è la larghezza di un inviluppo gaussiano nello spazio delle scale. Non è il raggio dello stencil della PDE, né il raggio di una singola sfocatura. È l’ampiezza del profilo di guadagno applicato alle bande wavelet discrete.

Lo schema di ponderazione è stabile indipendentemente dal livello di zoom : i raggi sono presi nello spazio dell’immagine raw (a piena risoluzione). Quando si visualizzano in anteprima immagini ridotte nella camera oscura, le frequenze più alte vengono troncate dal ridimensionamento, quindi applichiamo la decomposizione wavelet solo a partire dalla frequenza più alta disponibile e le ripesiamo in base al loro raggio equivalente a piena risoluzione. Ciò consente un’anteprima ridotta abbastanza accurata, anche se può nascondere artefatti di rumore che compaiono ai livelli più alti.

Le bande il cui raggio di sfocatura equivalente $\sigma_{G_s}$ si trova vicino a $r_c / z$ ricevono il guadagno maggiore, mentre le bande più lontane vengono progressivamente attenuate. Un $r_w$ piccolo dà una selezione ristretta di scale; un $r_w$ grande dà una risposta più ampia e più piatta sulle bande vicine.

Quindi la coppia (radius_center, radius) va letta come un centro e una larghezza nello spazio delle scale. Il modulo non definisce direttamente una gaussiana nella frequenza di Fourier; definisce invece un inviluppo gaussiano sulle bande multiscala disponibili.¹

Aggiornamento generico multiscala della PDE

Siano:

• $H_s$ la banda di dettaglio memorizzata all’indice di banda $s$,
• $G_s$ la ricostruzione a bassa frequenza corrente usata durante la risoluzione della banda $H_s$ nella sintesi.

A una data banda $s$, il codice costruisce quattro risposte di diffusione:

$$ \begin{align} D_{1,s} &= p_1 \, K_{1,s}(a_1, G_s) * G_s,\\ D_{2,s} &= p_2 \, K_{2,s}(a_2, H_s) * G_s, \\ D_{3,s} &= p_3 \, K_{3,s}(a_3, G_s) * H_s, \\ D_{4,s} &= p_4 \, K_{4,s}(a_4, H_s) * H_s, \end{align} $$

dove:

• $K_{1,s}$ e $K_{2,s}$ sono stencil di diffusione anisotropa 3x3 applicati alla ricostruzione a bassa frequenza corrente $G_s$,
• $K_{3,s}$ e $K_{4,s}$ sono gli stencil analoghi applicati alla banda di dettaglio $H_s$,
• da $a_1$ ad $a_4$ sono i 4 coefficienti di smorzamento dell’anisotropia definiti dall’utente (vedi sopra),
• da $p_1$ a $p_4$ sono i 4 coefficienti di aggiornamento (trasporto) della PDE definiti dall’utente (vedi sotto),
• la convoluzione è per canale sul sottoreticolo sparso à-trous di passo $2^s$.¹

Tutti i kernel laplaciani $K$ sono applicati con la stessa distanza di passo della sfocatura B-spline usata per produrre la scala wavelet $s$ su cui vengono applicati, cioè $2^s$. Dato che il kernel di sfocatura è 5×5 e i kernel laplaciani sono 3×3, ciò significa che il laplaciano copre un quarto della superficie della sfocatura.

Tutti questi $K$ sono uno sforzo per rendere il diffuse or sharpen un framework PDE multiscala generico, poiché incroceranno le informazioni tra $G_s$ e $H_s$ :

I gradienti valutati su $G_s$ hanno maggiore probabilità di puntare verso dettagli legittimi dell’immagine, e possono essere usati in un contesto di nitidezza per ignorare il rumore. I gradienti valutati su $H_s$ sono più sensibili al rumore e possono essere usati in un contesto diffusivo per attenuare il rumore. Separiamo il livello su cui inferiamo la struttura dell’immagine (direzione del gradiente) dal livello su cui calcoliamo l’aggiornamento della PDE (laplaciano) e quindi la diffusione.

In mancanza di un termine migliore, queste risposte $K_i$ sono legate a parametri della GUI chiamati «ordini», dal primo al quarto :

• speed, che è il coefficiente di aggiornamento della PDE,
• anisotropy, che è il coefficiente di smorzamento $a$ del tensore di anisotropia visto sopra.

Abbiamo perciò 4 parametri speed, uno per ciascun ordine :

$$ (p_1, p_2, p_3, p_4) = (\texttt{first}, \texttt{second}, \texttt{third}, \texttt{fourth}), $$

E analogamente, 4 parametri anisotropy : $$ (a_1, a_2, a_3, a_4) = (\texttt{first}, \texttt{second}, \texttt{third}, \texttt{fourth}) $$

Il significato di tutto ciò può essere tradotto in termini profani come segue :

1. diffondiamo la struttura nella direzione della struttura ($D_1$, primo ordine),
2. diffondiamo la struttura nella direzione della texture ($D_2$, secondo ordine),
3. diffondiamo la texture nella direzione della struttura ($D_3$, terzo ordine),
4. diffondiamo la texture nella direzione della texture ($D_4$, quarto ordine).

Qualsiasi coefficiente $p$ impostato a 0 annulla la diffusione, qualsiasi coefficiente $a$ impostato a 0 annulla l’anisotropia.

L’aggiornamento della PDE alla scala $s$ è perciò :

$$ U_s = G_s + H_s + \frac{\kappa \, w_s}{\nu_s} \sum_{i=1}^{4} D_{i,s}, $$

dove :

• $\kappa$ è il fattore di discretizzazione, cioè $\frac14$ per le differenze finite centrate,
• $\nu_s$ è il parametro di regolarizzazione che vedremo nella prossima sezione,
• $w_s$ è la ponderazione di scala definita nella sezione precedente.

E la risintesi finale è semplicemente :

$$ u’ = \sum_{s=0}^{n} U_s $$

Quindi, se riassumiamo l’intero algoritmo, dati $u$ l’immagine iniziale, $u’$ quella finale, $n$ il numero finale di scale :

$$ \begin{align*} s &\in [0, n] \\ G_{-1} &= u \\ G_s &= h_{2^s} * G_{s-1}, \\ H_s &= G_{s-1} - G_s, s\in[0,n] \\ w_s &= \exp\left( - \frac{(z\sigma_{G,s} - r_c)^2}{r_w^2} \right) \\ D_{1,s} &= p_1 \, K_{1,s}(a_1, G_s) * G_s \\ D_{2,s} &= p_2 \, K_{2,s}(a_2, H_s) * G_s \\ D_{3,s} &= p_3 \, K_{3,s}(a_3, G_s) * H_s \\ D_{4,s} &= p_4 \, K_{4,s}(a_4, H_s) * H_s \\ u’ &= \sum_{s=0}^{n} \left[ G_s + H_s + \frac14 \frac{w_s}{\nu_s} \sum_{i=1}^{4} D_{i,s}\right] \end{align*} $$

Alcune note:

• $n$ non è un parametro utente ma è determinato in funzione del $\sigma$ finale target richiesto dai parametri di raggio dell’utente. Come effetto collaterale, questo si adatta al livello di zoom.
• La sfocatura gaussiana è essa stessa la soluzione isotropa 2D dell’equazione del calore : è già diffusione.
• Aumentare il raggio di sfocatura equivale a lasciare che la diffusione prosegua per un tempo più lungo : il segnale si diffonde più lontano.
• Per $s > 0$, $H_s$ diventa in realtà una differenza di gaussiane. A un certo coefficiente di correzione di scala, la differenza di gaussiane è un’approssimazione del laplaciano di una gaussiana, che a sua volta è una stima del laplaciano a un coefficiente di scala $\sigma$.
• Applicare di nuovo un laplaciano (eventualmente anisotropo) su $H_s$ equivale a una derivata parziale del quarto ordine (bi-laplaciano).

Struttura del fattore di correzione

Il termine che controlla il trasporto della PDE è

$$ \frac{\kappa \, w_s}{\nu_s}, $$

In questa sezione definiremo il fattore $\nu_s$.

Definire filtri di frequenza ben comportati per le fotografie

Le fotografie sono riproduzioni digitali di un’immagine latente attraverso un apparato di esposizione (apertura del diaframma, sensibilità ISO del sensore, tempo di posa) e un apparato di discretizzazione (o campionamento spaziale) (matrice di filtri colore, griglia di pixel). Questi sono artefatti della tecnologia usata per catturare l’immagine latente e non riguardano l’immagine effettiva.

Sfortunatamente, le euristiche di cattura riguardanti esposizione e campionamento influenzano il modo in cui elaboriamo l’immagine digitale. Ho mostrato sui dischi nella sezione stencil di diffusione 3×3 come i bordi diagonali si comportino diversamente dai bordi allineati alla griglia (verticali/orizzontali), anche se ho scelto i kernel più invarianti alla rotazione : la rotazione del contenuto dell’immagine (rispetto alla griglia di pixel) cambierà il modo in cui i gradienti discreti vengono valutati. L’implicazione concreta qui è : ruotare l’immagine prima o dopo il diffuse or sharpen non produrrà lo stesso risultato.

Ma non finisce qui : l’esposizione cambia anche la varianza del segnale. Dato un segnale bianco $X$, la sua varianza locale $V_1$ su una finestra di campionamento $\mathcal{N}$ si esprime :

$$ \begin{align} V_{1, \mathcal{N}} &= \frac{1}{|\mathcal{N}|} \sum_{i\in\mathcal{N}} (\bar{X} - X_i)^2 \\ & = \frac{1}{|\mathcal{N}|} \sum_{i\in\mathcal{N}} \left(\left(\sum_{i\in\mathcal{N}} X_i \right) - X_i\right)^2 \end{align} $$

Se, invece di catturare $X$, sovraesponessimo la stessa immagine di un fattore $l$, allora la varianza di $lX$ diventerebbe $V_2$:

$$ \begin{align} V_{2, \mathcal{N}} &= \frac{1}{|\mathcal{N}|} \sum_{i\in\mathcal{N}} (l\bar{X} - lX_i)^2 \\ & = \frac{l^2}{|\mathcal{N}|} \sum_{i\in\mathcal{N}} (\bar{X} - X_i)^2 \\ & = l^2 \, V_{1, \mathcal{N}} \end{align} $$

Quindi la varianza del segnale aumenta con il quadrato del fattore di esposizione. Questo ci interessa qui per due ragioni, che possono essere riassunte fin d’ora così : tutto ciò che facciamo qui è cambiare la varianza dell’immagine scala per scala.

Innanzitutto, il nostro passo di sfocatura B-spline è una media locale pesata, e $H_s$ valutato al pixel di coordinate $(x, y)$ può in realtà essere scritto :

$$ \begin{align} H_s(x, y) &= X(x, y) - G_s(x,y)\\ &= X(x, y) - h_s * X(x,y)\\ &= -\left(\frac{1}{16^2} \sum_{i = -2}^{+2} \sum_{j = -2}^{+2} h_{i + 2} \, h_{j + 2} \, X(x + i, y + j) \right) + X(x, y) \end{align} $$

con $h_{i, j}$ i coefficienti del kernel B-spline 2D a 5 tap, e $X$ il segnale sfocato della scala precedente (o l’immagine iniziale per il primo passo). Possiamo mostrare in modo analogo che $H_s$ dipende linearmente dal fattore di esposizione. L’equazione sopra mostra come $H_s$ possa essere visto come una modulazione attorno a una media locale : la magnitudine di quella modulazione non è indipendente dalla magnitudine del segnale. Ciò significa che qualsiasi smussamento di $H_s$ (eseguito come processo diffusivo) avrà un peso diverso e un impatto diverso sui dettagli a seconda che l’immagine sia sovraesposta o sottoesposta, anche se il contenuto è lo stesso.

Detto altrimenti, smussare (o al contrario, aumentare la nitidezza) e poi sottoesporre, oppure sottoesporre e poi smussare non avranno lo stesso effetto sui dettagli, anche se la magnitudine complessiva finale (media) del segnale sarà la stessa. Non è ciò che ci aspettiamo da un filtro d’immagine ben comportato : la rappresentazione dei dati del contenuto non dovrebbe influenzare il modo in cui elaboriamo il contenuto stesso. Nel contesto diffusivo, non è così dannoso, ma nel contesto di nitidezza, i dettagli nelle ombre vengono davvero sovra-nitidizzati rispetto ai dettagli nelle alte luci, senza alcun tipo di normalizzazione.

In secondo luogo, la sfocatura B-spline (o la sua migliore approssimazione gaussiana) cambierà anche la varianza del segnale. Se esprimiamo il segnale discreto $X$ come una modulazione locale attorno alla media globale $\mu$, otteniamo $X_n = \mu + \epsilon_n$. Allora, la sfocatura B-spline applicata a $X$ diventa :

$$ \begin{align} G_{0}[n] &= \sum_{k} h_{k} (\mu + \epsilon_{n-k}) \\ &= \sum_{k} h_{k} \, \mu + \sum_{k} h_{k} \, \epsilon_{n-k} \\ &= \mu \sum_{k} h_{k} + \sum_{k} h_{k} \, \epsilon_{n-k} \end{align} $$

Poiché i coefficienti del kernel $h_{k}$ sono normalizzati e $\mu$ è, per definizione, costante sulla finestra di lunghezza $k$, otteniamo $\mu \sum_{k} h_{k} = \mu$, il che significa che la sfocatura non cambia il valore medio. La varianza è quindi espressa :

$$ \begin{align} \operatorname{Var}(G_0) &= \frac{1}{\mathcal{N}} \sum_{n \in \mathcal{N}} \left(\mu - \left(\mu + \sum_{k} h_{k} \, \epsilon_{n-k}\right) \right)^2 \\ &= \frac{1}{\mathcal{N}} \sum_{n \in \mathcal{N}} \left(\sum_{k} h_{k} \, \epsilon_{n-k} \right)^2 \\ &= \operatorname{Var}\left(\sum_{k} h_{k} \, \epsilon_{n-k}\right) \end{align} $$

Da qui, possiamo mostrare che, per segnali bianchi e non correlati, $\operatorname{Var}(G_0) = \operatorname{Var}(X) \sum_k h_k^2$, e più in generale :

$$ \operatorname{Var}(G_s) = \operatorname{Var}(G_{s-1}) \sum_k h_k^2 $$

dove $\sum_k h_k^2 = (35 / 128)^2$ per il filtro B-spline cardinale 2D a 5 tap.

E anche qui abbiamo un problema : lo stesso oggetto campionato a una certa risoluzione, o a 4 volte quella risoluzione, apparirebbe alla stessa «frequenza» un passo di decomposizione wavelet dopo, il che renderebbe la sua varianza di superficie inferiore di $(35 / 128)^2$. Ma questo è puramente un artefatto di campionamento. La varianza dell’oggetto stesso può essere concettualizzata, al di fuori dell’immagine, in un contesto continuo, come le modulazioni locali di colore attorno al suo colore di superficie medio. Sebbene questa varianza idealizzata non possa essere recuperata da alcun apparato di imaging, il modo in cui trattiamo il segnale dovrebbe almeno essere stabile in varianza, così da poter postulare che la varianza dell’immagine rappresenti la varianza idealizzata dell’oggetto a meno di un fattore di scala costante.

Questo potrebbe sembrare una preoccupazione filosofica finché non ci imbattiamo in un problema pratico di qualsiasi software di imaging : cosa succede quando visualizziamo l’anteprima dell’effetto ingrandendo/riducendo lo zoom ? Come scaliamo l’effetto affinché l’anteprima ridimensionata sia ancora fedele al risultato a piena risoluzione ?

Quindi, tutte queste discrepanze di campionamento devono essere normalizzate per ottenere un filtro d’immagine che tenti di manipolare il contenuto indipendentemente dalla sua rappresentazione dei dati.

Definire una metrica di regolarizzazione

Abbiamo visto sopra come la varianza del segnale sia una metrica rilevante per ciò che stiamo facendo qui : possiamo tracciarla attraverso i passi di sfocatura, collegarla alla magnitudine del segnale, ed essa rappresenta la modulazione del segnale attorno al valore medio.

Sfortunatamente, non abbiamo accesso a una metrica di varianza una volta entrati nello schema di decomposizione wavelet. Tuttavia, abbiamo visto sopra che $H_s$ era piuttosto vicino, concettualmente, al termine $(\bar{X} - X_i)$ della varianza :

• invece di una media aritmetica, usiamo una media pesata con i coefficienti B-spline,
• invece di una media globale, ne usiamo una locale,
• la natura radiale della B-spline la rende più invariante alla rotazione di qualsiasi media per patch quadrata.

Quindi useremo l’energia della banda $H_s$, valutata alle stesse coordinate di pixel dello stencil laplaciano, definita come :

$$ Q_s = \sum_{q \in \mathcal{N}_{3\times 3}} H_s(q)^2. $$

Per un segnale bianco, non correlato e a variazione lenta, $\overline{Q_s} = Q_s / |\mathcal{N}_{3\times 3}|$ diventa prossimo alla varianza per scala e per patch.

La regolarizzazione è pensata per il problema della nitidezza, che è mal definito : in questo contesto, aumentiamo l’energia di ciascun livello $H_s$ e ci serve un parametro per tenerla a bada a un certo punto. Questa è una procedura comune nei problemi inversi come la riduzione del rumore e la deconvoluzione, per i quali la Variazione Totale è stata usata come schema di regolarizzazione per parecchio tempo.

Il modello di regolarizzazione che useremo è :

$$ \nu_s = \tau + \lambda \, \dfrac{1}{9} \sum_{q \in \mathcal{N}_{3\times 3}} \left(\frac{H_s(q)}{L_s(q)} \right)^2. $$

con i parametri utente :

$$ \lambda = 10^{\texttt{regularization}} - 1, \qquad \tau = 10^{\texttt{variance_threshold}}. $$

Abbiamo mostrato sopra come la varianza del segnale vari con il quadrato del fattore di scala dell’esposizione, e come $H_s$ vari linearmente con il fattore di scala dell’esposizione. $G_s$ porta la stessa dipendenza lineare attraverso la sua proprietà di essere una media locale pesata.

Quindi, il rapporto $H_s / G_s$ è invariante rispetto all’esposizione. Per identificazione, $L_s(q) = G_s(q)$ nell’equazione di regolarizzazione, quindi usiamo l’energia di banda invariante rispetto all’esposizione :

$$ Q_s’ = \sum_{q \in \mathcal{N}_{3\times 3}} \left(\frac{H_s(q)}{G_s(q)} \right)^2 $$

e la sua media locale :

$$ \overline{Q_s’} = \frac{1}{9} \sum_{q \in \mathcal{N}_{3\times 3}} \left(\frac{H_s(q)}{G_s(q)} \right)^2 $$

Normalizzare scala e copertura spaziale

Lo stencil laplaciano 3×3 si espande di un passo $2^s$ per ogni scala $s$, proprio come fa il kernel B-spline : questa è la base dello schema «à-trous». Lo spazio fisico coperto da quel kernel aumenta con le scale wavelet.

Possiamo mostrare che un filtro gaussiano 2D di parametro $\sigma$ cambia la varianza del segnale originale come :

$$ \operatorname{Var}(g_\sigma * X) = \dfrac{\operatorname{Var}(X)}{4 \pi \sigma^2} $$

Dato che $4 \pi \sigma^2$ è l’area effettiva del disco coperta dal filtro gaussiano, ciò corrisponde a un’intuizione molto semplice : il filtro gaussiano diffonde la modulazione originale (espressa come varianza) su una superficie più ampia. Questa è la diffusione in poche parole. Da ciò deriviamo :

$$ \operatorname{Var}(X) \propto \sigma^2 \operatorname{Var}(g_\sigma * X) $$

Tra 2 passi di sfocatura gaussiana, il parametro di varianza del filtro gaussiano equivalente varia di :

$$ \begin{align} \Delta \sigma_s^2 &= \sigma_s^2 - \sigma_{s-1}^2 \\ &= \sigma_B^2 \, 4^s \end{align} $$

E da una scala all’altra, questo raggio cresce di $\Delta \sigma = \sigma_B 2^s \approx 2^s$, ricordando che $\sigma_B \approx 1.05\dots$.

Quindi la valutazione del laplaciano si espande spazialmente allo stesso ritmo (stesso passo a ogni scala) dei filtri B-spline o gaussiani : il laplaciano segue implicitamente la diffusione della varianza attraverso le scale, e qui non è necessaria alcuna normalizzazione di scala aggiuntiva.

L’implementazione iniziale del diffuse or sharpen aveva un potenziamento $\sigma^2_s$ applicato al parametro di regolarizzazione $\lambda$ : in pratica, ciò impedisce alle scale grossolane di avere qualsiasi impatto visibile, e le limita a essere il supporto per il lavoro alle alte frequenze.

È stato fatto un altro tentativo per normalizzare $\lambda$ tenendo conto del fatto che $Q_s$ ha un’energia di banda decrescente al crescere di $s$, quindi una metrica di energia invariante di scala portava a :

$$ E_s = \frac{4}{(\Delta\sigma_s^2)^2}Q_s $$

Anche combinato con il fattore $\sigma^2_s$ sopra, ciò dava un peso di gran lunga eccessivo alle scale grossolane, rendendo difficile aumentare la nitidezza alle alte frequenze mentre quelle basse stavano già sovrastimando pesantemente (dando un aspetto impastato). Entrambi questi tentativi sono stati abbandonati.

Risultato

After

Before

L’immagine qui sopra sembrerà eccessiva ai più, ma non è questo il punto: ho avuto la mia dose di operatori di nitidezza difettosi che funzionavano abbastanza bene finché non spingevi l’intensità oltre il 2%. Filtri d’immagine di questo tipo non hanno alcun interesse. Impari molto di più sugli algoritmi osservando come falliscono che come riescono nel loro punto ottimale.

Questa immagine è la trappola definitiva per ogni algoritmo di nitidezza:

• il primo piano è più vicino e meno velato dello sfondo, quindi implora di essere sovra-nitidizzato,
• il primo piano ravvicinato è anche molto più scuro, quindi di nuovo è facile spingere la nitidezza verso il cielo lì,
• abbiamo una cresta montuosa molto contrastata che implora di produrre aloni intorno ai bordi,
• il disco del sole in un cielo nuvoloso verrà nitidizzato dalla maggior parte degli algoritmi con bordi scuri,
• la quantità di rimozione della foschia necessaria farebbe esplodere il rumore (questa foto è stata scattata nel 2017 con una Nikon D90 a 200 ISO, siamo lontani dai sensori attuali).

Quindi, qui eseguiamo una riduzione del rumore e una rimozione della sfocatura congiunte a grandi raggi. La rimozione della sfocatura usa una contro-diffusione isotropica (che descrive al meglio la foschia atmosferica), sia sulle alte frequenze sia sulle basse frequenze. La riduzione del rumore usa la diffusione lungo gli isofoti sulle alte frequenze, seguendo il gradiente campionato nelle alte frequenze. Tutto questo è stato fatto senza alcuna maschera, in una singola istanza di diffuse or sharpen.

Nel complesso, non vediamo alcun sovrascatto dei bordi né aloni. Il primo piano scuro è stato per lo più ignorato, come dovrebbe essere, ma abbiamo recuperato i dettagli giù nella valle. Non abbiamo ottenuto alcuno spostamento di colore né aberrazioni cromatiche. Non pretenderò che questa immagine sia priva di artefatti, perché, guardando da vicino, otteniamo striature di colore e nuovi dettagli che potrebbero essere una ricostruzione reale o pure allucinazioni del modello diffusivo. Ma resta il fatto che questi artefatti, se presenti, appaiono abbastanza organici da passare inosservati se non abbiamo l’immagine originale.

Questo convalida la pertinenza dello schema diffusivo multiscala, insieme alla strategia di regolarizzazione.

After

Before

After

Before

Questa è stata scattata con una Sony Alpha 7R Mk 2 ma con un Konika Hexanon 40 mm f/1.8. È un obiettivo pancake degli anni ‘70, e quindi la foto è composta da 45 Mpx di splendida sfocatura. I ritratti sono meno indulgenti dei paesaggi perché la pelle deve rimanere sana quando si recupera la nitidezza. Qui abbiamo anche tutto il lens flare e i raggi solari che possono degenerare rapidamente.

Il dettaglio è particolarmente impressionante qui perché recuperiamo le sopracciglia e le ciglia senza sovra-nitidizzare i capelli in controluce. Nel complesso, non c’è ringing, fringing o aloni. Purtroppo, la qualità della pelle si è degradata ed è il massimo che possiamo raggiungere senza mascherare la pelle per escluderla. Da qui in poi, si tratta solo di riportare selettivamente il contrasto globale con una curva tonale:

After

Before

Cosa è specifico di Ansel

I riferimenti scientifici spiegano i mattoni costitutivi:

• diffusione anisotropica e inpainting per trasferimento di calore,²
• scale-space e diffusione discreta,³⁴
• laplaciani isotropici a 9 punti,⁵⁶
• approssimazione gaussiana con B-spline.⁷

Ciò che è specifico di Ansel è il modo in cui sono assemblati:

• una piramide B-spline à-trous a piena risoluzione,
• quattro operatori di diffusione parametrizzati indipendentemente, suddivisi tra basse e alte frequenze,
• un regolarizzatore dell’energia di banda HF normalizzato inoltre in base all’energia LF locale,
• selezione della scala consapevole dello zoom nell’anteprima della camera oscura,
• matematica identica su CPU e OpenCL.¹

Quindi il modulo va inteso come una sintesi ingegneristica di diverse idee numeriche, non come l’implementazione letterale di un singolo articolo.

Prospettive

Lo schema wavelet multiscala di PDE di diffusione anisotropica con regolarizzazione produce risultati fotografici utilizzabili, al di là della semplice prova di concetto. Permette di sfruttare la nitidezza e la riduzione del rumore congiunte, insieme alla consueta diffusione orientata. Permette anche di ringiovanire, tramite software, vecchi obiettivi ritenuti inadatti alla fotografia digitale ad alta definizione. Il modulo stesso fornisce un ambiente generico per le PDE che può essere usato in molti modi diversi.

Il problema, però, è che la natura delle impostazioni è radicata in matematica di livello (almeno) universitario ed è criptica per la maggior parte dei fotografi. Spiegare l’effetto dei parametri è difficile senza addentrarsi in ciò che significano matematicamente. Cercare di rinominarli in base alla loro funzione anziché alla loro natura è destinato a fallire perché la loro funzione dipende da come vengono combinati tra loro, e dal fatto che vengano usati nell’intervallo di valori positivi o negativi.

La configurazione della diffusione è piuttosto semplice e non richiede regolarizzazione. Non ha bisogno di tutti e 4 gli ordini insieme, ma bastano solo le impostazioni del primo e del terzo ordine. Nella configurazione isotropica, è del tutto equivalente a una sfocatura gaussiana, che sarà meno costosa da calcolare (perché non iterativa).

La configurazione della nitidezza, insieme alla riduzione del rumore congiunta, è più complicata. L’unico modo per renderla più intuitiva è addestrarla come un algoritmo di apprendimento automatico:

• scattare coppie di immagini pulite e sfocate/rumorose/velate della stessa scena (sfocatura da movimento, sfocatura da sfocatura dell’obiettivo, obiettivi morbidi),
• eseguire una scansione dei parametri a forza bruta del modulo diffuse or sharpen che esegue la ricostruzione delle immagini sporche, e registrare la norma $L_2$ dell’errore tra le immagini di riferimento pulite e le ricostruzioni tentate,
• classificare manualmente le immagini in categorie binarie (sfocata, rumorosa, velata) o per intensità (il rumore può essere misurato come PSNR, RMS, …; probabilmente anche il raggio dovrà essere presente),
• il problema di apprendimento automatico diventa: nello spazio 14D dei parametri di ingresso di diffuse or sharpen, quali sono le 4 direzioni principali (associate a riduzione del rumore, rimozione della sfocatura, rimozione della foschia, raggio) che minimizzano l’errore $E = ||\text{clean} - \text{reconstructed}||_2$? Cerchiamo gli autovettori di quelle direzioni principali.
• risolverlo tramite PLS pesato (partial least-squares): $Y = X’ B + c$ per $n$ coppie di immagini pulite/ricostruite, dove:
  • $Y$ è il vettore $4 × n$ delle categorie per ogni campione,
  • $X’$ è il vettore $14 × n$ standardizzato (standardizzato componente per componente: $X_j’ = \frac{X_r - \mu}{\sigma}$) dei parametri di diffuse or sharpen,
  • $B$ è la matrice $14 × 4$ che mappa i nostri 14 criptici parametri di D or S su 4 parametri intuitivi (è l’incognita qui),
  • $c$ è il residuo (scalare o vettore a seconda di ciò che si adatta),
  • la norma $L_2$ dell’errore è usata come ponderazione della PLS (probabilmente iniettata in un esponenziale),
• una volta nota la matrice $B$ (mappa 14D -> 4D), invertirla per ottenere il modello 4D -> 14 D,
• aggiungere una modalità GUI alternativa che espone i 4 parametri intuitivi, e uno strato GUI <-> parametri che converte quei 4 nei 14 argomenti di ingresso di D or S (ovvero: scrivere il prodotto matriciale).

Qualsiasi altro tentativo di “semplificare” diffuse or sharpen sarà solo un’insensata operazione di rietichettatura che offusca il significato reale dei parametri e impedisce a chiunque abbia le adeguate basi matematiche di comprenderlo. Sarebbe un peccato fare in modo che le poche persone in grado di comprenderlo vengano di fatto dissuase persino dal provarci. Al momento, D or S è difficile da comprendere, ma almeno può essere spiegato. Rietichettare i controlli non lo renderà più facile da comprendere, ma aggiungerà solo un ulteriore strato di traduzione semantica tra matematica e GUI, molto probabilmente comunque impreciso e fuorviante, il che lo renderà solo più impegnativo dal punto di vista cognitivo da spiegare e da afferrare.