De wiskunde van diffuse or sharpen

Dit artikel documenteert de wiskunde die momenteel is geïmplementeerd in Ansels diffuse or sharpen-module, zoals te vinden in src/iop/diffuse.c, src/common/bspline.h, data/kernels/diffuse.cl en data/kernels/bspline.cl. Het is geen gebruikshandleiding. Het is een reconstructie van het numerieke model op basis van de broncode, waarbij de wetenschappelijke beweringen worden herleid tot de referenties die in de codecommentaren worden aangehaald.¹

Samenvatting

De module werkt op een multischaaldecompositie van het beeld op volledige resolutie volgens het à-trous-schema, opgebouwd uit kardinale B-spline-vervagingen. Op elke schaal berekent Ansel vier anisotrope tweede-orde-operatoren op de laagfrequente en hoogfrequente banden afzonderlijk, reguleert hun werking met een op schaal genormaliseerde lokale hoogfrequente bandenergie, beperkt de update optioneel tot een binair inpainting-masker en reconstrueert het beeld van grof naar fijn.² De discrete ruimtelijke operatoren zijn 3x3-stencils met centrale verschillen, bemonsterd op het ijle à-trous-subrooster van de huidige schaal, zodat de grotere fysieke drager voortkomt uit de stap $2^s$ van de waveletladder in plaats van uit een grotere PDE-kern.³¹

Motivatie

Het doel van de diffuse-module was om te simuleren hoe aquarelachtig pigment uit een beeld naar de randen ervan ontsnapt, zoals dit:

Het bleek dat het simpelweg omkeren van het teken van de update van de partiële differentiaalvergelijking voor diffusie het beeld daadwerkelijk kon verscherpen. Ik had jarenlang aan blinde deconvolutie gewerkt toen ik met een diffusiemodel begon, en het lukte me nooit om zelfs de state-of-the-art-methode buiten de mooie voorbeeldgevallen (lage ISO, lage verzadiging, uniforme vervaging) te laten werken. Maar dit had natuurlijk een regularisatieparameter nodig om divergentie te voorkomen.

Zo werd de diffuse or sharpen-module geboren, als een generiek raamwerk voor het simuleren van allerlei diffusieve of contradiffusieve verschijnselen.

Continu model

In de kern implementeert de module een familie van anisotrope diffusievergelijkingen van de vorm

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

waarbij $u$ het beeld is en $\mathbf{A}$ een symmetrische positieve diffusietensor die het gladstrijken stuurt, ofwel langs de isofoten, langs de gradiënt, of isotroop.²⁴

De oorspronkelijke inspiratie in de code is het inpainting-model van Qin et al., die beeldstructuur en textuurherstel koppelen via anisotrope warmteoverdracht.² Ansel behoudt die geest, maar past de PDE toe in een waveletdomein, afzonderlijk op laagfrequente en hoogfrequente banden, en stelt vier door de gebruiker bestuurbare transportcoëfficiënten beschikbaar.

Discrete gradiënten en anisotropie

Voor elke pixel en elk kanaal extraheert de module een 3x3-omgeving en evalueert gecentreerde eerste afgeleiden

$$ \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} $$

Dit is de standaard gradiënt met centrale verschillen die wordt gebruikt in discrete schaalruimte- en diffusiemethoden.³

De lokale gradiëntoriëntatie is dan:

$$ \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} $$

met de gebruikelijke terugval op $(1,0)$ wanneer de gradiëntmagnitude verdwijnt.¹ Merk op dat we het dure berekenen van de hoek $\theta$ als de $\arctan2$-functie van de $(x, y)$-componenten zullen vermijden, aangezien deze in het volgende nooit rechtstreeks zal worden gebruikt.

De anisotropiesterkte wordt omgezet van de gebruikersparameter $a$ naar een positieve coëfficiënt

$$ \alpha = a^2, $$

terwijl het teken van $a$ de modus selecteert:

• $a = 0$: isotrope diffusie,
• $a > 0$: diffusie uitgelijnd met de isofootrichting,
• $a < 0$: diffusie uitgelijnd met de gradiëntrichting.¹

De dempingsterm is dan

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

wat overeenkomt met de anisotrope diffusie van Qin et al. voor inpainting.²

Voor isofoot-uitgelijnde diffusie is de tensor, geschreven in de lokale beeldbasis:

$$ \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}. $$

Voor gradiënt-uitgelijnde diffusie gebruikt Ansel de omgekeerde vorm:

$$ \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}. $$

De 3x3-diffusiestencils

Afhankelijk van de gebruikersparameters zullen we voor elke pixel het $\mathbf{K(A)}$-stencil evalueren met behulp van een van de bovenstaande $\mathbf{A}$-tensoren.

Gegeven de symmetrische tensor $\mathbf{A}$:

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

bouwen we de discrete, geroteerde, anisotrope Laplace-kern $\nabla \cdot (\mathbf{A}\nabla u)$ als het 3x3-stencil:

$$ \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}, $$

op de tekenconventie na die in de broncode wordt vermeld voor de niet-diagonale termen.²⁴¹

In het isotrope geval verwordt dit tot een constante:

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

wat de klassieke isotrope Laplace-operator in Oono-Puri-stijl is. Het belang ervan ligt in het rotatiegedrag: vergeleken met de 5-punts-Laplace-operator is de hoekfout kleiner, wat van belang is wanneer diffusie de beeldassen niet mag bevoordelen.⁵⁶ Deze hoekfout werd geëvalueerd ten opzichte van vergelijkbare kernen in Rotation-invariant Laplacian for 2D grids .

Om het effect van richting in diffusie te begrijpen, laten we beginnen met een ruizige schijf en die met Ansel diffunderen:

Een paar dingen zijn hier het vermelden waard:

• Alleen de isotrope vervaging levert een achromatische schijf op, de andere behouden wat grootschalige chromaruis, omdat de chromaruis lokale gradiëntvariaties veroorzaakt.
• De anisotrope diffusie is perfect nabij de noord/oost/zuid/west-posities van de cirkel, dat wil zeggen wanneer de gradiënthoek perfect is uitgelijnd op het pixelraster (0° of 90°). Daartussenin zien we afwijkingen in de manier waarop met richting wordt omgegaan, als gevolg van de numerieke beperkingen van een geroteerd vierkant 3×3-stencil.

De à-trous-B-spline-piramide

De multischaalanalyse is opgebouwd uit het scheidbare 5-taps kardinale B-spline-filter

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

Dit filter is een compacte benadering van een gaussische kern, en daarom verschijnt het van nature in spline-gebaseerde schaalruimtemethoden.⁷ We hergebruiken hier het raamwerk dat door Johannes Hanika werd geïntroduceerd voor de contrast equalizer-module.⁸

De huidige implementatie gebruikt de historische ongedecimeerde à-trous-ladder:

1. De eerste laagdoorlaat op volledige resolutie is

   $$ G_0 = h_0 * u, $$

2. De fijnste detailband is

   $$ H_0 = u - G_0. $$

3. Grovere niveaus behouden dezelfde beeldresolutie en vergroten alleen de vervagingsstap met $2^s$:

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

   wat betekent dat bij $s=1$ de $h_{2}$-taps om de andere pixel worden geëvalueerd, bij $s=2$ de $h_4$-taps om de 4 pixels worden geëvalueerd, enzovoort…

4. De band die op schaal $s$ wordt opgeslagen is het verschil tussen twee opeenvolgende laagdoorlaten:

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

Zo leeft elke band op het oorspronkelijke beeldraster (geen decimatie), en lossen we de PDE fijn-naar-grof parallel op, door de oplossing schaal voor schaal in de uitvoerbuffer te accumuleren totdat we het uiteindelijke residu toevoegen.¹ Dit schema vermijdt het opstapelen van de fijne-schaalcorrectie bovenop de grove schalen, zoals gebruikelijk is bij gaussische piramides of multigrid-oplossers, en heeft empirisch een betere stabiliteit aangetoond bij het inverse probleem van scherpteherstel. Dit staat gelijk aan het afzonderlijk beheren van de energie van elke band, op een manier die lijkt op HiFi-audio-equalizers, maar dan binnen een 2D-ruimtelijk raamwerk.

Kernschaal, bandschaal en de GUI-schaalenvelop

Om de eigenschappen van het waveletschema te bestuderen, zullen we gebruikmaken van de eigenschap van de kardinale B-spline dat deze een benadering van een gaussisch filter is. De gaussische parameter $\sigma$ bepaalt hoe vervagingskernen samenstellen:

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

dus gaussische kernvarianties $\sigma^2$ tellen op bij convolutie.³⁷ Daarentegen hangt de variantie van een gefilterd signaal $X$:

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

af van het spectrum van $X$ en is over het algemeen niet gelijk aan $\sigma^2$. Om beide varianties te onderscheiden, verwijst $\sigma$ in de rest van dit artikel alleen naar de parameter van de gaussische kern, nooit naar de vierkantswortel van de signaalvariantie.

Om dubbelzinnigheid te vermijden, helpt het ook om de laagdoorlaatniveaus te scheiden van de niveaubanden:

• $G_0 = h_1 * u$ is de eerste laagdoorlaat op volledige resolutie over het invoerbeeld $u$,
• $G_s = h_{2^s} * G_{s-1}$ voor $s \ge 1$,
• $H_0 = u - G_0$,
• $H_s = G_{s-1} - G_s$ voor $s \ge 1$.

Dus:

• de laagdoorlaatindex $s$ verwijst naar het vervagingsniveau $G_s$,
• de bandindex $s$ verwijst naar de detailband $H_s$ die tussen $G_{s-1}$ en $G_s$ ligt,
• elk niveau blijft bemonsterd op het oorspronkelijke beeldraster.¹

Het B-spline-filter benadert het best de equivalente gaussische kern met parameter:

$$ \sigma_B \approx 1.055365. $$

Omdat gaussische kernvarianties optellen bij convolutie, volgt de effectieve vervagingsstraal van $G_s$ rechtstreeks uit de à-trous-analyseladder. Daarom is de cumulatieve laagdoorlaatreeks

$$ 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 $$

met

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

en, in het algemeen,

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

Equivalent,

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

De code gebruikt $\sigma_{G_s}$ om te bepalen hoeveel schalen er nodig zijn om overeen te komen met de door de gebruiker gevraagde radius-parameter en om elke band rond de door de gebruiker geselecteerde centrale straal te wegen:

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

waarbij $z$ het zoomniveau van de donkere kamer is (vergeleken met raw op volledige resolutie), $r_c$ is radius_center in de GUI, en $r_w$ is radius.¹ Door deze in te stellen, definiëren gebruikers het verval van een soort gediscretiseerd “bandbreedte”-filter gecentreerd op een willekeurige frequentie. In die zin is $r_w$ de breedte van een gaussische envelop in de schaalruimte. Het is niet de straal van het PDE-stencil, noch de straal van een enkele vervaging. Het is de spreiding van het versterkingsprofiel dat op de discrete waveletbanden wordt toegepast.

Het wegingsschema is stabiel ongeacht het zoomniveau: stralen worden genomen in de raw-beeldruimte (volledige resolutie). Bij het bekijken van in de donkere kamer verkleinde beelden worden de hoogste frequenties door het verkleinen afgekapt, dus passen we de waveletdecompositie alleen toe vanaf de hoogste beschikbare frequentie en herwegen we deze op basis van hun equivalente straal op volledige resolutie. Dit maakt een tamelijk nauwkeurige verkleinde voorvertoning mogelijk, ook al kan het ruizige artefacten verbergen die op de hoogste niveaus verschijnen.

Banden waarvan de equivalente vervagingsstraal $\sigma_{G_s}$ dicht bij $r_c / z$ ligt, krijgen de grootste versterking, terwijl banden die verder weg liggen progressief worden verzwakt. Een kleine $r_w$ geeft een smalle selectie van schalen; een grote $r_w$ geeft een bredere en vlakkere respons over naburige banden.

Het paar (radius_center, radius) moet dus worden gelezen als een centrum en een breedte in de schaalruimte. De module definieert niet rechtstreeks een Gauss in de Fourier-frequentie; in plaats daarvan definieert ze een Gaussische omhullende over de beschikbare multischaalbanden.¹

Generieke multischaal-PDE-update

Laat:

• $H_s$ de opgeslagen detailband bij bandindex $s$ zijn,
• $G_s$ de huidige laagfrequente reconstructie zijn die wordt gebruikt bij het oplossen van band $H_s$ tijdens de synthese.

Bij een gegeven band $s$ bouwt de code vier diffusieresponsen op:

$$ \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} $$

waarbij:

• $K_{1,s}$ en $K_{2,s}$ 3x3 anisotrope diffusiestencils zijn die worden toegepast op de huidige laagfrequente reconstructie $G_s$,
• $K_{3,s}$ en $K_{4,s}$ de analoge stencils zijn die worden toegepast op de detailband $H_s$,
• $a_1$ tot $a_4$ 4 door de gebruiker gedefinieerde anisotropie-dempingscoëfficiënten zijn (zie boven),
• $p_1$ tot $p_4$ 4 door de gebruiker gedefinieerde PDE-update-coëfficiënten (transport) zijn (zie onder),
• de convolutie kanaalsgewijs is op het ijle à-trous-subrooster met stapgrootte $2^s$.¹

Alle Laplaciaanse $K$-kernels worden toegepast met dezelfde stapafstand als de B-spline-vervaging die is gebruikt om de waveletschaal $s$ te produceren waarover zij worden toegepast, dat wil zeggen $2^s$. Aangezien de vervagingskernel 5×5 is en de Laplaciaanse kernels 3×3 zijn, betekent dit dat de Laplaciaan een kwart van het vervagingsoppervlak beslaat.

Al die $K$ zijn een poging om diffuse or sharpen tot een generiek multischaal-PDE-raamwerk te maken, aangezien zij informatie zullen kruisen tussen $G_s$ en $H_s$ :

De gradiënten die op $G_s$ worden geëvalueerd, wijzen waarschijnlijker naar legitieme beelddetails en kunnen in een verscherpingscontext worden gebruikt om ruis te negeren. De gradiënten die op $H_s$ worden geëvalueerd, zijn gevoeliger voor ruis en kunnen in een diffusieve context worden gebruikt om ruis te verzachten. We scheiden de laag waarop we de beeldstructuur afleiden (gradiëntrichting) van de laag waar we de PDE-update berekenen (Laplaciaan) en dus de diffusie.

Bij gebrek aan een betere term zijn die $K_i$-responsen gekoppeld aan GUI-parameters die „orders" (ordes) worden genoemd, van de eerste tot de vierde :

• speed, wat de PDE-update-coëfficiënt is,
• anisotropy, wat de $a$-dempingscoëfficiënt uit de anisotropietensor hierboven is.

We hebben dus 4 speed-parameters, één voor elke orde :

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

En op vergelijkbare wijze 4 anisotropy-parameters : $$ (a_1, a_2, a_3, a_4) = (\texttt{first}, \texttt{second}, \texttt{third}, \texttt{fourth}) $$

De betekenis van dit alles kan in lekentermen als volgt worden vertaald :

1. we diffunderen structuur in de richting van structuur ($D_1$, eerste orde),
2. we diffunderen structuur in de richting van textuur ($D_2$, tweede orde),
3. we diffunderen textuur in de richting van structuur ($D_3$, derde orde),
4. we diffunderen textuur in de richting van textuur ($D_4$, vierde orde).

Elke $p$-coëfficiënt die op 0 is ingesteld, heft diffusie op, elke $a$-coëfficiënt die op 0 is ingesteld, heft anisotropie op.

De PDE-update op schaal $s$ is daarom :

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

waarbij :

• $\kappa$ de discretisatiefactor is, dat wil zeggen $\frac14$ voor eindige gecentreerde differenties,
• $\nu_s$ de regularisatieparameter is die we in de volgende sectie zullen zien,
• $w_s$ de schaalweging is die in de vorige sectie is gedefinieerd.

En de uiteindelijke hersynthese is eenvoudigweg :

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

Dus, als we het hele algoritme samenvatten, gegeven $u$ het initiële beeld, $u’$ het uiteindelijke, $n$ het uiteindelijke aantal schalen :

$$ \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*} $$

Enkele opmerkingen:

• $n$ is geen gebruikersparameter maar wordt bepaald met betrekking tot de beoogde uiteindelijke $\sigma$ die door de gebruikersstraalparameters wordt gevraagd. Dit past zich als bijproduct aan het zoomniveau aan.
• De Gaussische vervaging is zelf de 2D-isotrope oplossing van de warmtevergelijking : het is al diffusie.
• Het vergroten van de vervagingsstraal is equivalent aan het langer laten lopen van de diffusie : het signaal verspreidt zich verder.
• Voor $s > 0$ wordt $H_s$ feitelijk een verschil van Gaussianen. Op een bepaalde schalingscorrectiecoëfficiënt is het verschil van Gaussianen een benadering van de Laplaciaan van een Gauss, die zelf een schatting is van de Laplaciaan op een schalingscoëfficiënt $\sigma$.
• Het opnieuw toepassen van een (mogelijk anisotrope) Laplaciaan over $H_s$ is equivalent aan een vierde-orde partiële afgeleide (bi-Laplaciaan).

Structuur van de correctiefactor

De term die het PDE-transport regelt is

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

In deze sectie zullen we de $\nu_s$-factor definiëren.

Goed gedragen frequentiefilters voor foto’s definiëren

Foto’s zijn digitale reproducties van een latent beeld via een belichtingsapparaat (diafragma-opening, ISO-gevoeligheid van de sensor, sluitertijd) en een discretisatie- (of ruimtelijke bemonsterings-)apparaat (kleurenfilterarray, pixelraster). Dit zijn artefacten van de technologie die wordt gebruikt om het latente beeld vast te leggen en betreffen niet het eigenlijke beeld.

Helaas beïnvloeden de opnameheuristieken met betrekking tot belichting en bemonstering hoe we het digitale beeld verwerken. Ik heb in de sectie 3×3 diffusiestencils met schijven laten zien hoe diagonale randen zich anders gedragen dan rasteruitgelijnde randen (verticaal/horizontaal), ook al koos ik de meest rotatie-invariante kernels : rotatie van de beeldinhoud (ten opzichte van het pixelraster) verandert hoe discrete gradiënten worden geëvalueerd. De concrete implicatie hier is : het beeld roteren vóór of na diffuse or sharpen zal niet hetzelfde resultaat opleveren.

Maar het houdt daar niet op : belichting verandert ook de signaalvariantie. Gegeven een wit signaal $X$, wordt zijn lokale variantie $V_1$ over een bemonsteringsvenster $\mathcal{N}$ uitgedrukt als :

$$ \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} $$

Als we, in plaats van $X$ vast te leggen, hetzelfde beeld met een factor $l$ zouden overbelichten, dan zou de variantie van $lX$ $V_2$ worden:

$$ \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} $$

De variantie van het signaal neemt dus toe met het kwadraat van de belichtingsfactor. Dit is voor ons hier van belang om twee redenen, die nu al als volgt kunnen worden samengevat : alles wat we hier doen is het veranderen van de schaalsgewijze variantie van het beeld.

Ten eerste is onze B-spline-vervagingsstap een lokaal gewogen gemiddelde, en $H_s$ geëvalueerd bij de pixel met coördinaten $(x, y)$ kan feitelijk worden geschreven als :

$$ \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} $$

met $h_{i, j}$ de coëfficiënten van de 2D B-spline-kernel met 5 taps, en $X$ het vervaagde signaal van de vorige schaal (of het initiële beeld voor de eerste stap). We kunnen op vergelijkbare wijze aantonen dat $H_s$ lineair afhankelijk is van de belichtingsfactor. De bovenstaande vergelijking laat zien hoe $H_s$ kan worden gezien als een modulatie rond een lokaal gemiddelde : de amplitude van die modulatie is niet onafhankelijk van de amplitude van het signaal. Dit betekent dat elke gladstrijking van $H_s$ (uitgevoerd als een diffusief proces) een ander gewicht en een andere impact op details zal hebben, afhankelijk van of het beeld over- of onderbelicht is, ook al is de inhoud dezelfde.

Anders gezegd, gladstrijken (of omgekeerd, verscherpen) en dan onderbelichten, of onderbelichten en dan gladstrijken zal niet hetzelfde effect op details hebben, ook al zal de uiteindelijke totale amplitude (gemiddelde) van het signaal hetzelfde zijn. Dit is niet wat we verwachten van een goed gedragen beeldfilter : de datarepresentatie van de inhoud zou niet mogen beïnvloeden hoe we de inhoud zelf verwerken. In de diffusieve context is dat niet zo schadelijk, maar in de verscherpingscontext worden details in schaduwen echt overmatig verscherpt vergeleken met details in hooglichten, zonder enige vorm van normalisatie.

Ten tweede zal de B-spline-vervaging (of zijn beste Gaussische benadering) ook de variantie van het signaal veranderen. Als we het discrete signaal $X$ uitdrukken als een lokale modulatie rond het globale gemiddelde $\mu$, krijgen we $X_n = \mu + \epsilon_n$. Dan wordt de B-spline-vervaging toegepast op $X$ :

$$ \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} $$

Omdat de coëfficiënten van de $h_{k}$-kernel genormaliseerd zijn en $\mu$ per definitie constant is over het venster van lengte $k$, krijgen we $\mu \sum_{k} h_{k} = \mu$, wat betekent dat de vervaging de gemiddelde waarde niet verandert. De variantie wordt dan uitgedrukt als :

$$ \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} $$

Vanaf daar kunnen we aantonen dat, voor witte, ongecorreleerde signalen, $\operatorname{Var}(G_0) = \operatorname{Var}(X) \sum_k h_k^2$, en meer algemeen :

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

waarbij $\sum_k h_k^2 = (35 / 128)^2$ voor het 2D kardinale B-spline-filter met 5 taps.

En ook hier hebben we een probleem : hetzelfde object, bemonsterd op een bepaalde resolutie, of op 4 keer die resolutie, zou één stap van waveletdecompositie later op dezelfde „frequentie" verschijnen, wat zijn oppervlaktevariantie $(35 / 128)^2$ lager zou maken. Maar dit is puur een bemonsteringsartefact. De variantie van het object zelf kan worden geconceptualiseerd, buiten het beeld, in een continue context, als de lokale kleurmodulaties rond zijn gemiddelde oppervlaktekleur. Hoewel deze geïdealiseerde variantie door geen enkel beeldvormingsapparaat kan worden hersteld, zou de manier waarop we het signaal behandelen ten minste stabiel moeten zijn in variantie, zodat we kunnen stellen dat de beeldvariantie de geïdealiseerde variantie van het object representeert op een constante schalingsfactor na.

Dit lijkt misschien een filosofische bekommernis totdat we op een praktisch probleem van elke beeldsoftware stuiten : wat gebeurt er wanneer we het effect in-/uitgezoomd bekijken ? Hoe schalen we het effect zodat de verkleinde voorbeeldweergave nog steeds getrouw is aan het resultaat op volle resolutie ?

Dus al deze bemonsteringsdiscrepanties moeten worden genormaliseerd om een beeldfilter te bereiken dat inhoud probeert te manipuleren ongeacht zijn datarepresentatie.

Een regularisatiemetriek definiëren

We hebben hierboven gezien hoe de signaalvariantie een relevante metriek is voor wat we hier doen : we kunnen ze door de vervagingsstappen volgen, ze koppelen aan de signaalamplitude, en ze representeert de signaalmodulatie rond de gemiddelde waarde.

Helaas hebben we geen toegang tot een variantiemetriek zodra we het waveletdecompositieschema binnengaan. We hebben echter hierboven gezien dat $H_s$ conceptueel vrij dicht bij de $(\bar{X} - X_i)$-term van de variantie lag :

• in plaats van een rekenkundig gemiddelde gebruiken we een gewogen gemiddelde met B-spline-coëfficiënten,
• in plaats van een globaal gemiddelde gebruiken we een lokaal gemiddelde,
• het radiale karakter van de B-spline maakt ze rotatie-invarianter dan enig vierkant patchgewijs gemiddelde.

Dus zullen we de energie van de $H_s$-band gebruiken, geëvalueerd bij dezelfde pixelcoördinaten als de Laplaciaanse stencil, gedefinieerd als :

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

Voor een langzaam variërend, wit, ongecorreleerd signaal wordt $\overline{Q_s} = Q_s / |\mathcal{N}_{3\times 3}|$ dicht bij de schaalsgewijze en patchgewijze variantie.

De regularisatie is bedoeld voor het verscherpingsprobleem, dat slecht gedefinieerd is : in deze context verhogen we de energie van elke laag $H_s$ en hebben we een parameter nodig om ze op een gegeven moment in toom te houden. Dit is een gebruikelijke procedure bij inverse problemen zoals ruisonderdrukking en ontwazing, waarvoor Total Variation al geruime tijd als regularisatieschema is gebruikt.

Het regularisatiemodel dat we zullen gebruiken is :

$$ \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. $$

met de gebruikersparameters :

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

We hebben hierboven aangetoond hoe de signaalvariantie varieert met het kwadraat van de belichtingsschaling, en hoe $H_s$ lineair varieert met de belichtingsschaling. $G_s$ draagt dezelfde lineaire afhankelijkheid via zijn eigenschap een lokaal gewogen gemiddelde te zijn.

Dus de verhouding $H_s / G_s$ is belichtingsinvariant. Door identificatie is $L_s(q) = G_s(q)$ in de regularisatievergelijking, dus gebruiken we de belichtingsinvariante bandenergie :

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

en zijn lokale gemiddelde :

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

Schaal en ruimtelijke dekking normaliseren

De 3×3 Laplaciaanse stencil breidt zich uit met een stapgrootte van $2^s$ voor elke schaal $s$, net zoals de B-spline-kernel doet : dit is de basis van het „à-trous"-schema. De fysieke ruimte die door die kernel wordt bestreken, neemt toe met de waveletschalen.

We kunnen aantonen dat een 2D Gaussisch filter met parameter $\sigma$ de oorspronkelijke signaalvariantie verandert als :

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

Aangezien $4 \pi \sigma^2$ het effectieve schijfoppervlak is dat door het Gaussische filter wordt bestreken, komt dit overeen met een zeer eenvoudige intuïtie : het Gaussische filter spreidt de oorspronkelijke modulatie (uitgedrukt als variantie) over een groter oppervlak. Dat is diffusie in een notendop. Daaruit leiden we af :

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

Tussen 2 Gaussische vervagingsstappen varieert de variantieparameter van het equivalente Gaussische filter met :

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

En van de ene schaal naar de andere groeit deze straal met $\Delta \sigma = \sigma_B 2^s \approx 2^s$, in herinnering brengend dat $\sigma_B \approx 1.05\dots$.

Dus de evaluatie van de Laplaciaan breidt zich ruimtelijk uit met hetzelfde tempo (dezelfde stapgrootte bij elke schaal) als de B-spline- of Gaussische filters : de Laplaciaan volgt impliciet de variantiespreiding over de schalen, en er is hier geen aanvullende schaalnormalisatie nodig.

De initiële implementatie van diffuse or sharpen had een $\sigma^2_s$-boost toegepast op de $\lambda$-regularisatieparameter : in de praktijk voorkomt dat dat grove schalen enige zichtbare impact hebben, en beperkt het ze tot het zijn van de ondersteuning voor hoogfrequent werk.

Er werd nog een poging gedaan om $\lambda$ te normaliseren, rekening houdend met het feit dat $Q_s$ een afnemende bandenergie heeft naarmate $s$ toeneemt, dus een schaal-invariante energiemetriek leidde tot :

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

Zelfs gecombineerd met de $\sigma^2_s$-factor hierboven, gaf dat veel te veel gewicht aan de grove schalen, waardoor het moeilijk werd om bij hoge frequenties te verscherpen terwijl de lage al ernstig overschoten (wat een uitgesmeerde look gaf). Beide pogingen zijn opgegeven.

Resultaat

After

Before

De afbeelding hierboven zal er voor de meesten overdreven uitzien, maar dat is niet waar het om gaat: ik heb mijn portie kapotte verscherpingsoperatoren gehad die goed genoeg werkten zolang je de sterkte niet boven de 2% opdreef. Zulke beeldfilters zijn totaal oninteressant. Je leert veel meer over algoritmes door te kijken naar hoe ze falen dan naar hoe ze slagen in hun ideale werkgebied.

Deze afbeelding is de ultieme valstrik voor elk verscherpingsalgoritme:

• de voorgrond is dichterbij en minder heiig dan de achtergrond, dus die smeekt erom te worden oververscherpt,
• de nabije voorgrond is ook veel donkerder, dus opnieuw is het makkelijk om de scherpte daar de lucht in te jagen,
• we hebben een zeer contrastrijke bergkam die smeekt om halo’s rond randen te produceren,
• de zonneschijf in een bewolkte lucht zal door de meeste algoritmes met donkere randen worden verscherpt,
• de hoeveelheid ontwaseming die nodig is zou de ruis doen exploderen (deze is genomen in 2017 met een Nikon D90 op 200 ISO, we zijn ver verwijderd van de huidige sensoren).

Dus hier voeren we een gezamenlijke ruisonderdrukking en ontwazing uit bij grote radii. De ontwazing gebruikt isotrope tegendiffusie (die atmosferische heiigheid het best beschrijft), op zowel hoogfrequentie als laagfrequentie. De ruisonderdrukking gebruikt isofoot-diffusie op hoogfrequentie, volgens de in hoogfrequentie bemonsterde gradiënt. Dit alles werd zonder enige maskering gedaan, in één enkele instantie van diffuse or sharpen.

Al met al zien we geen randoverschot en geen halo’s. De donkere voorgrond werd grotendeels genegeerd, zoals het hoort, maar we kregen de details in het dal terug. We kregen geen kleurverschuiving of chromatische aberraties. Ik zal niet beweren dat deze afbeelding artefactvrij is, want als je van dichtbij kijkt, krijgen we kleurstrepen en nieuwe details die daadwerkelijke reconstructie of pure hallucinaties van het diffusiemodel kunnen zijn. Maar het punt blijft dat deze artefacten, als ze er al zijn, er organisch genoeg uitzien om onopgemerkt te blijven als we niet over de originele afbeelding beschikken.

Dit bevestigt de relevantie van het meerschalige diffusieschema, samen met de regularisatiestrategie.

After

Before

After

Before

Deze is genomen met een Sony Alpha 7R Mk 2 maar met een Konika Hexanon 40 mm f/1.8. Dit is een pancake-lens uit de jaren 1970, en dus is de afbeelding 45 Mpx aan prachtige onscherpte. Portretten zijn minder vergevingsgezind dan landschappen omdat de huid gezond moet blijven wanneer je de scherpte terugbrengt. Hier hebben we ook alle lensflare en de zonnestralen die snel kunnen ontaarden.

Het detail is hier bijzonder indrukwekkend omdat we de wenkbrauwen en wimpers terugkrijgen zonder het tegenlichthaar te oververscherpen. Al met al is er geen ringing, fringing of haloing. Helaas is de huidkwaliteit verslechterd en dit is zo ver als we kunnen gaan zonder de huid te maskeren om die uit te sluiten. Vanaf daar gaat het er alleen nog om het globale contrast selectief terug te brengen met een tooncurve:

After

Before

Wat specifiek is voor Ansel

De wetenschappelijke referenties leggen de bouwstenen uit:

• anisotrope diffusie en warmteoverdracht-inpainting,²
• schaalruimte en discrete diffusie,³⁴
• isotrope 9-punts Laplacianen,⁵⁶
• B-spline Gaussische benadering.⁷

Wat specifiek is voor Ansel is de manier waarop ze worden samengesteld:

• een a-trous B-spline-piramide op volledige resolutie,
• vier onafhankelijk geparametriseerde diffusieoperatoren, verdeeld over lage en hoge frequenties,
• een HF-bandenergie-regularisator die aanvullend genormaliseerd wordt door de lokale LF-energie,
• zoombewuste schaalselectie in de voorbeeldweergave van de donkere kamer,
• identieke CPU- en OpenCL-wiskunde.¹

De module moet dus worden opgevat als een technische synthese van meerdere numerieke ideeën, niet als een letterlijke implementatie van één enkel artikel.

Perspectieven

Het meerschalige waveletschema van anisotrope diffusie-PDE met regularisatie levert bruikbare fotografische resultaten op, verder dan alleen het proof of concept. Het maakt het mogelijk om gezamenlijke verscherping en ruisonderdrukking te benutten, samen met reguliere georiënteerde diffusie. Het maakt het ook mogelijk om, via software, oude lenzen te verjongen die ongeschikt werden geacht voor hoge-definitie digitale fotografie. De module zelf biedt een generieke PDE-speeltuin die op veel verschillende manieren kan worden gebruikt.

Het probleem is echter dat de aard van de instellingen geworteld is in (minstens) wiskunde op bachelorniveau en cryptisch is voor de meeste fotografen. Het effect van de parameters uitleggen is moeilijk zonder in te duiken in wat ze wiskundig betekenen. Ze proberen te hernoemen naar hun functie in plaats van hun aard is gedoemd te mislukken, omdat hun functie afhangt van hoe ze met elkaar worden gecombineerd, en of ze in het positieve of negatieve waardenbereik worden gebruikt.

De diffusieopzet is vrij rechttoe rechtaan en vereist geen regularisatie. Ze heeft niet alle 4 de ordes tegelijk nodig, maar alleen de instellingen van de eerste en derde orde volstaan. In de isotrope instelling is ze volledig equivalent aan een Gaussische vervaging, die goedkoper te berekenen is (omdat niet-iteratief).

De verscherpingsopzet, samen met gezamenlijke ruisonderdrukking, is ingewikkelder. De enige manier om die gebruiksvriendelijker te maken is door haar te trainen als een machine-learning-algoritme:

• schiet paren van schone en onscherpe/ruizige/heiige foto’s van dezelfde scène (bewegingsonscherpte, lensdefocusonscherpte, zachte lenzen),
• voer een brute-force parameter-sweep uit van de diffuse or sharpen-module die reconstructie van de vuile afbeeldingen uitvoert, en registreer de $L_2$-norm van de fout tussen de schone referentieafbeeldingen en de gepoogde reconstructies,
• classificeer de afbeeldingen handmatig in binaire categorieën (onscherp, ruizig, heiig) of naar intensiteit (ruis kan worden gemeten als PSNR, RMS, …; radius zal er waarschijnlijk ook bij moeten zijn),
• het machine-learning-probleem wordt: in de 14D-ruimte van diffuse or sharpen-invoerparameters, wat zijn de 4 hoofdrichtingen (geassocieerd met ruisonderdrukking, ontwazing, ontwaseming, radius) die de fout $E = ||\text{clean} - \text{reconstructed}||_2$ minimaliseren? We zoeken naar de eigenvectoren van die hoofdrichtingen.
• los dat op met gewogen PLS (partiële kleinste kwadraten): $Y = X’ B + c$ voor $n$ paren schone/gereconstrueerde afbeeldingen, waarbij:
  • $Y$ de $4 × n$-vector van categorieën voor elk monster is,
  • $X’$ de $14 × n$ gestandaardiseerde vector is (component-gewijs gestandaardiseerd: $X_j’ = \frac{X_r - \mu}{\sigma}$) van diffuse or sharpen-parameters,
  • $B$ de $14 × 4$-matrix is die onze 14 cryptische D or S-parameters afbeeldt op 4 gebruiksvriendelijke parameters (dit is hier de onbekende),
  • $c$ het residu is (scalair of vector, afhankelijk van wat past),
  • de $L_2$-norm van de fout wordt gebruikt als de PLS-weging (waarschijnlijk geïnjecteerd in een exponentiële functie),
• zodra de matrix $B$ bekend is (afbeelding 14D -> 4D), inverteer die om het 4D -> 14D-model te verkrijgen,
• voeg een alternatieve GUI-modus toe die de 4 gebruiksvriendelijke parameters blootstelt, en een GUI <-> parameters-laag die die 4 omzet naar de 14 invoerargumenten van D or S (betekent: schrijf het matrixproduct).

Elke andere poging om diffuse or sharpen te “vereenvoudigen” zal slechts een dwaze her-etiketteringsklus zijn die de werkelijke betekenis van de parameters verhult en verhindert dat iemand met de juiste wiskundige achtergrond het begrijpt. Het zou jammer zijn om ervoor te zorgen dat de enige weinigen die het kunnen begrijpen juist ontmoedigd zouden worden om het zelfs maar te proberen. Op dit moment is D or S moeilijk te begrijpen, maar het kan tenminste worden uitgelegd. Bedieningselementen her-etiketteren zal het niet makkelijker te begrijpen maken, maar zal alleen een extra laag van semantische vertaling tussen wiskunde en GUI toevoegen, hoogstwaarschijnlijk toch onnauwkeurig en misleidend, wat het alleen maar cognitief veeleisender zal maken om uit te leggen en te bevatten.