Die Mathematik von Diffundieren oder Schärfen

Dieser Artikel dokumentiert die Mathematik, die derzeit in Ansels Modul diffuse or sharpen implementiert ist, wie sie in src/iop/diffuse.c, src/common/bspline.h, data/kernels/diffuse.cl und data/kernels/bspline.cl zu finden ist. Es ist keine Anwendungsanleitung. Es ist eine Rekonstruktion des numerischen Modells aus dem Quellcode, wobei die wissenschaftlichen Aussagen auf die in den Code-Kommentaren zitierten Referenzen zurückgeführt werden.¹

Kurzfassung

Das Modul operiert auf einer À-trous-Mehrskalen-Zerlegung des Bildes in voller Auflösung, die aus kardinalen B-Spline-Unschärfen aufgebaut ist. Auf jeder Skala berechnet Ansel vier anisotrope Operatoren zweiter Ordnung getrennt auf dem Niederfrequenz- und dem Hochfrequenzband, reguliert deren Wirkung durch eine skalennormierte lokale Hochfrequenzband-Energie, beschränkt die Aktualisierung optional auf eine binäre Inpainting-Maske und rekonstruiert das Bild von grob zu fein.² Die diskreten räumlichen Operatoren sind 3x3-Stencils zentraler Differenzen, die auf dem dünn besetzten À-trous-Untergitter der aktuellen Skala abgetastet werden, sodass die größere physikalische Reichweite aus der Schrittweite $2^s$ der Wavelet-Leiter statt aus einem größeren PDE-Kern stammt.³¹

Motivation

Das Ziel des Diffuse-Moduls war es, zu simulieren, wie aquarellartiges Pigment aus einem Bild in seine Ränder entweicht, etwa so :

Es stellte sich heraus, dass die bloße Umkehrung des Vorzeichens der Aktualisierung der partiellen Differentialgleichung der Diffusion das Bild tatsächlich schärfen konnte. Ich hatte jahrelang an blinder Dekonvolution gearbeitet, als ich anfing, an einem Diffusionsmodell zu arbeiten, und es nie geschafft, selbst die dem Stand der Technik entsprechende Methode außerhalb der schönen Vorzeigefälle (niedrige ISO, geringe Sättigung, gleichmäßige Unschärfe) zum Funktionieren zu bringen. Aber dies erforderte natürlich einen Regularisierungsparameter, um ein Divergieren zu vermeiden.

So entstand das Modul diffuse or sharpen, als generisches Framework für die Simulation aller Arten von diffusiven oder gegendiffusiven Phänomenen.

Kontinuierliches Modell

Im Kern implementiert das Modul eine Familie anisotroper Diffusionsgleichungen der Form

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

wobei $u$ das Bild ist und $\mathbf{A}$ ein symmetrischer positiver Diffusionstensor, der die Glättung entweder entlang der Isophoten, entlang des Gradienten oder isotrop lenkt.²⁴

Die ursprüngliche Inspiration im Code ist das Inpainting-Modell von Qin et al., die Bildstruktur und Texturrekonstruktion durch anisotropen Wärmetransfer koppeln.² Ansel behält diesen Geist bei, wendet die PDE aber in einem Wavelet-Bereich an, getrennt auf Niederfrequenz- und Hochfrequenzbändern, und stellt vier benutzergesteuerte Transportkoeffizienten bereit.

Diskrete Gradienten und Anisotropie

Für jedes Pixel und jeden Kanal extrahiert das Modul eine 3x3-Nachbarschaft und wertet zentrierte erste Ableitungen aus

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

Dies ist der Standardgradient zentraler Differenzen, der in diskreten Skalenraum- und Diffusionsmethoden verwendet wird.³

Die lokale Gradientenorientierung ist dann :

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

mit dem üblichen Rückgriff auf $(1,0)$, wenn der Gradientenbetrag verschwindet.¹ Man beachte, dass wir die Berechnung des aufwendigen Winkels $\theta$ als $\arctan2$-Funktion der Komponenten $(x, y)$ vermeiden werden, da er im Folgenden nie direkt verwendet wird.

Die Anisotropiestärke wird vom Benutzerparameter $a$ in einen positiven Koeffizienten umgewandelt

$$ \alpha = a^2, $$

während das Vorzeichen von $a$ den Modus auswählt:

• $a = 0$: isotrope Diffusion,
• $a > 0$: Diffusion ausgerichtet entlang der Isophotenrichtung,
• $a < 0$: Diffusion ausgerichtet entlang der Gradientenrichtung.¹

Der Dämpfungsterm ist dann

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

was der anisotropen Diffusion von Qin et al. für Inpainting entspricht.²

Für isophotenausgerichtete Diffusion lautet der in der lokalen Bildbasis geschriebene Tensor :

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

Für gradientenausgerichtete Diffusion verwendet Ansel die invertierte Form :

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

Die 3x3-Diffusions-Stencils

Abhängig von den Benutzerparametern werten wir für jedes Pixel das Stencil $\mathbf{K(A)}$ unter Verwendung eines der obigen Tensoren $\mathbf{A}$ aus.

Gegeben den symmetrischen Tensor $\mathbf{A}$ :

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

bauen wir den diskreten, rotierten, anisotropen Laplace-Kern $\nabla \cdot (\mathbf{A}\nabla u)$ als das 3x3-Stencil auf :

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

bis auf die im Quellcode erwähnte Vorzeichenkonvention für die Außerdiagonalterme.²⁴¹

Im isotropen Fall vereinfacht es sich zu einer Konstanten :

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

was der klassische isotrope Laplace-Operator im Oono-Puri-Stil ist. Sein Nutzen liegt im Rotationsverhalten: im Vergleich zum 5-Punkt-Laplace-Operator wird der Winkelfehler reduziert, was wichtig ist, wenn die Diffusion die Bildachsen nicht bevorzugen darf.⁵⁶ Dieser Winkelfehler wurde gegen ähnliche Kerne in Rotation-invariant Laplacian for 2D grids  ausgewertet.

Um die Wirkung der Richtung bei der Diffusion zu verstehen, beginnen wir mit einer verrauschten Scheibe und diffundieren sie mit Ansel :

Ein paar Dinge sind hier erwähnenswert :

• Nur die isotrope Unschärfe ergibt eine achromatische Scheibe, die anderen behalten ein gewisses großskaliges Chroma-Rauschen, weil das Chroma-Rauschen lokale Gradientenvariationen erzeugt.
• Die anisotrope Diffusion ist nahe den Nord-/Ost-/Süd-/West-Positionen des Kreises perfekt, das heißt, wenn der Gradientenwinkel perfekt am Pixelgitter ausgerichtet ist (0° oder 90°). Dazwischen sehen wir Abweichungen in der Behandlung der Richtung, bedingt durch die numerischen Beschränkungen eines rotierten quadratischen 3×3-Stencils.

Die À-trous-B-Spline-Pyramide

Die Mehrskalenanalyse wird aus dem separierbaren 5-Tap-Kardinal-B-Spline-Filter aufgebaut

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

Dieser Filter ist eine kompakte Gauß-Approximation, weshalb er in spline-basierten Skalenraummethoden auf natürliche Weise auftritt.⁷ Wir verwenden hier das von Johannes Hanika für das Modul contrast equalizer eingeführte Framework wieder.⁸

Die aktuelle Implementierung verwendet die historische, nicht dezimierte À-trous-Leiter:

1. Der erste Tiefpass in voller Auflösung ist

   $$ G_0 = h_0 * u, $$

2. Das feinste Detailband ist

   $$ H_0 = u - G_0. $$

3. Gröbere Ebenen behalten dieselbe Bildauflösung und vergrößern nur die Unschärfe-Schrittweite um $2^s$:

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

   das heißt, bei $s=1$ werden die Taps $h_{2}$ an jedem nächsten Pixel ausgewertet, bei $s=2$ werden die Taps $h_4$ an jedem 4. Pixel ausgewertet, und so weiter…

4. Das auf Skala $s$ gespeicherte Band ist die Differenz zwischen zwei aufeinanderfolgenden Tiefpässen:

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

Somit lebt jedes Band auf dem ursprünglichen Bildgitter (keine Dezimierung), und wir lösen die PDE parallel von fein zu grob, indem wir die Lösung Skala für Skala in den Ausgabepuffer akkumulieren, bis wir das endgültige Residuum hinzufügen.¹ Dieses Schema vermeidet das Aufstapeln der Korrektur feiner Skalen auf die groben Skalen, wie es bei Gauß-Pyramiden oder Mehrgitter-Lösern üblich ist, und hat empirisch eine bessere Stabilität beim inversen Problem der Schärferekonstruktion gezeigt. Dies entspricht der getrennten Verwaltung der Energie jedes Bandes, in einer Weise, die HiFi-Audio-Equalizern ähnelt, aber innerhalb eines 2D-räumlichen Rahmens.

Kern-Skala, Band-Skala und die Skalenhüllkurve der GUI

Um die Eigenschaften des Wavelet-Schemas zu untersuchen, nutzen wir die Eigenschaft des kardinalen B-Splines, eine Approximation eines Gauß-Filters zu sein. Der Gauß-Parameter $\sigma$ steuert, wie sich Unschärfekerne zusammensetzen:

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

sodass sich die Varianzen $\sigma^2$ der Gauß-Kerne unter Faltung addieren.³⁷ Im Gegensatz dazu hängt die Varianz eines gefilterten Signals $X$ :

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

vom Spektrum von $X$ ab und ist im Allgemeinen nicht gleich $\sigma^2$. Um beide Varianzen zu unterscheiden, bezieht sich $\sigma$ im Rest dieses Artikels nur auf den Parameter des Gauß-Kerns, niemals auf die Quadratwurzel der Signalvarianz.

Um Mehrdeutigkeit zu vermeiden, hilft es außerdem, die Tiefpassebenen von den Bändern der Ebenen zu trennen:

• $G_0 = h_1 * u$ ist der erste Tiefpass in voller Auflösung über dem Eingabebild $u$,
• $G_s = h_{2^s} * G_{s-1}$ für $s \ge 1$,
• $H_0 = u - G_0$,
• $H_s = G_{s-1} - G_s$ für $s \ge 1$.

Also:

• der Tiefpassindex $s$ bezieht sich auf die Unschärfeebene $G_s$,
• der Bandindex $s$ bezieht sich auf das Detailband $H_s$, das zwischen $G_{s-1}$ und $G_s$ liegt,
• jede Ebene bleibt auf dem ursprünglichen Bildgitter abgetastet.¹

Der B-Spline-Filter approximiert am besten den äquivalenten Gauß-Kern mit dem Parameter :

$$ \sigma_B \approx 1.055365. $$

Da sich die Varianzen der Gauß-Kerne unter Faltung addieren, folgt der effektive Unschärferadius von $G_s$ direkt aus der À-trous-Analyseleiter. Daher ist die kumulative Tiefpassfolge

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

mit

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

und, allgemein,

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

Gleichwertig,

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

Der Code verwendet $\sigma_{G_s}$, um zu entscheiden, wie viele Skalen benötigt werden, um dem vom Benutzer angeforderten Parameter radius zu entsprechen, und um jedes Band um den vom Benutzer gewählten zentralen Radius zu gewichten:

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

wobei $z$ die Zoomstufe der Dunkelkammer ist (im Vergleich zum RAW in voller Auflösung), $r_c$ ist radius_center in der GUI und $r_w$ ist radius.¹ Durch das Festlegen dieser Werte definieren Benutzer den Abfall einer Art diskretisierten „Bandbreiten"-Filters, der auf einer beliebigen Frequenz zentriert ist. In diesem Sinne ist $r_w$ die Breite einer Gauß-Hüllkurve im Skalenraum. Es ist nicht der Radius des PDE-Stencils, noch der Radius einer einzelnen Unschärfe. Es ist die Streuung des Verstärkungsprofils, das auf die diskreten Wavelet-Bänder angewendet wird.

Das Gewichtungsschema ist unabhängig von der Zoomstufe stabil : die Radien werden im RAW-Bildraum (volle Auflösung) genommen. Bei der Vorschau von in der Dunkelkammer herunterskalierten Bildern werden die höchsten Frequenzen durch die Herunterskalierung abgeschnitten, sodass wir die Wavelet-Zerlegung nur ab der höchsten verfügbaren Frequenz anwenden und sie entsprechend ihrem äquivalenten Radius in voller Auflösung neu gewichten. Dies ermöglicht eine recht genaue herunterskalierte Vorschau, auch wenn sie verrauschte Artefakte verbergen kann, die auf den höchsten Ebenen auftreten.

Bänder, deren äquivalenter Unschärferadius $\sigma_{G_s}$ nahe bei $r_c / z$ liegt, erhalten die größte Verstärkung, während weiter entfernte Bänder fortschreitend abgeschwächt werden. Ein kleines $r_w$ ergibt eine schmale Auswahl von Skalen; ein großes $r_w$ ergibt eine breitere und flachere Antwort über benachbarte Bänder.

Das Paar (radius_center, radius) ist also als Zentrum und Breite im Skalenraum zu lesen. Das Modul definiert keinen Gauß direkt in der Fourier-Frequenz; stattdessen definiert es eine Gaußsche Einhüllende über die verfügbaren Multiskalen-Bänder.¹

Generisches Multiskalen-PDE-Update

Sei:

• $H_s$ das gespeicherte Detailband beim Bandindex $s$,
• $G_s$ die aktuelle Niederfrequenz-Rekonstruktion, die beim Lösen von Band $H_s$ während der Synthese verwendet wird.

Bei einem gegebenen Band $s$ erzeugt der Code vier Diffusionsantworten:

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

wobei:

• $K_{1,s}$ und $K_{2,s}$ 3x3-Stencils der anisotropen Diffusion sind, die auf die aktuelle Niederfrequenz-Rekonstruktion $G_s$ angewandt werden,
• $K_{3,s}$ und $K_{4,s}$ die analogen Stencils sind, die auf das Detailband $H_s$ angewandt werden,
• $a_1$ bis $a_4$ 4 benutzerdefinierte Anisotropie-Dämpfungskoeffizienten sind (siehe oben),
• $p_1$ bis $p_4$ 4 benutzerdefinierte PDE-Update-Koeffizienten (Transportkoeffizienten) sind (siehe unten),
• die Faltung kanalweise auf dem dünn besetzten à-trous-Untergitter mit Schrittweite $2^s$ erfolgt.¹

Alle Laplace-Kerne $K$ werden mit derselben Schrittweite (stride) angewandt wie die B-Spline-Unschärfe, die zur Erzeugung der Wavelet-Skala $s$ diente, über die sie angewandt werden, also $2^s$. Da der Unschärfekern 5×5 groß ist und die Laplace-Kerne 3×3, bedeutet das, dass der Laplace-Kern ein Viertel der Unschärfefläche abdeckt.

Alle diese $K$ sind ein Versuch, diffuse or sharpen zu einem generischen Multiskalen-PDE-Framework zu machen, denn sie tauschen Information zwischen $G_s$ und $H_s$ aus:

Die auf $G_s$ ausgewerteten Gradienten zeigen mit höherer Wahrscheinlichkeit in Richtung legitimer Bilddetails und können in einer Schärfungssituation genutzt werden, um Rauschen zu ignorieren. Die auf $H_s$ ausgewerteten Gradienten sind rauschempfindlicher und können in einer diffusiven Situation genutzt werden, um Rauschen zu glätten. Wir trennen die Ebene, auf der wir die Bildstruktur ableiten (Gradientenrichtung), von der Ebene, auf der wir das PDE-Update (Laplace) und damit die Diffusion berechnen.

Mangels eines besseren Begriffs sind diese $K_i$-Antworten mit GUI-Parametern verknüpft, die „orders" (Ordnungen) genannt werden, von der ersten bis zur vierten:

• speed, das der PDE-Update-Koeffizient ist,
• anisotropy, das der Dämpfungskoeffizient $a$ aus dem Anisotropie-Tensor weiter oben ist.

Wir haben daher 4 speed-Parameter, einen für jede Ordnung:

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

Und analog 4 anisotropy-Parameter: $$ (a_1, a_2, a_3, a_4) = (\texttt{first}, \texttt{second}, \texttt{third}, \texttt{fourth}) $$

Die Bedeutung all dessen lässt sich für Laien etwa so übersetzen:

1. wir diffundieren Struktur in Richtung der Struktur ($D_1$, erste Ordnung),
2. wir diffundieren Struktur in Richtung der Textur ($D_2$, zweite Ordnung),
3. wir diffundieren Textur in Richtung der Struktur ($D_3$, dritte Ordnung),
4. wir diffundieren Textur in Richtung der Textur ($D_4$, vierte Ordnung).

Jeder auf 0 gesetzte $p$-Koeffizient hebt die Diffusion auf, jeder auf 0 gesetzte $a$-Koeffizient hebt die Anisotropie auf.

Das PDE-Update bei Skala $s$ lautet daher:

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

wobei:

• $\kappa$ der Diskretisierungsfaktor ist, also $\frac14$ für zentrierte finite Differenzen,
• $\nu_s$ der Regularisierungsparameter ist, den wir im nächsten Abschnitt sehen werden,
• $w_s$ die Skalengewichtung ist, die im vorherigen Abschnitt definiert wurde.

Und die endgültige Resynthese lautet einfach:

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

Wenn wir also den gesamten Algorithmus zusammenfassen, gegeben $u$ das Ausgangsbild, $u’$ das endgültige, $n$ die endgültige Anzahl der Skalen:

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

Einige Anmerkungen:

• $n$ ist kein Benutzerparameter, sondern wird im Hinblick auf das vom Benutzer über die Radiusparameter angeforderte endgültige Ziel-$\sigma$ bestimmt. Dies passt sich als Nebeneffekt an den Zoomgrad an.
• Die Gaußsche Unschärfe ist selbst die isotrope 2D-Lösung der Wärmeleitungsgleichung: sie ist bereits Diffusion.
• Das Vergrößern des Unschärferadius ist äquivalent dazu, die Diffusion für eine längere Zeit laufen zu lassen: das Signal breitet sich weiter aus.
• Für $s > 0$ wird $H_s$ tatsächlich zu einer Differenz von Gaußfunktionen. Bei einem gewissen Skalierungskorrekturkoeffizienten ist die Differenz von Gaußfunktionen eine Näherung des Laplace einer Gaußfunktion, der selbst wiederum eine Schätzung des Laplace bis auf einen Skalierungskoeffizienten $\sigma$ ist.
• Das erneute Anwenden eines (möglicherweise anisotropen) Laplace über $H_s$ ist äquivalent zu einer partiellen Ableitung vierter Ordnung (Bi-Laplace).

Struktur des Korrekturfaktors

Der Term, der den PDE-Transport steuert, ist

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

In diesem Abschnitt werden wir den Faktor $\nu_s$ definieren.

Definition gutartiger Frequenzfilter für Fotografien

Fotografien sind digitale Reproduktionen eines latenten Bildes durch einen Belichtungsapparat (Blendenöffnung, ISO-Empfindlichkeit des Sensors, Verschlusszeit) und einen Diskretisierungsapparat (bzw. räumliche Abtastung) (Farbfilter-Array, Pixelraster). Das sind Artefakte der zur Erfassung des latenten Bildes verwendeten Technik und betreffen nicht das eigentliche Bild.

Leider beeinflussen die Erfassungsheuristiken bezüglich Belichtung und Abtastung, wie wir das digitale Bild verarbeiten. Ich habe im Abschnitt 3×3-Diffusions-Stencils an Scheiben gezeigt, wie sich diagonale Kanten anders verhalten als rasterausgerichtete Kanten (vertikal/horizontal), obwohl ich die rotationsinvariantesten Kerne gewählt habe: die Rotation des Bildinhalts (gegenüber dem Pixelraster) verändert, wie diskrete Gradienten ausgewertet werden. Die konkrete Folge hier ist: das Bild vor oder nach diffuse or sharpen zu drehen, ergibt nicht dasselbe Ergebnis.

Aber damit nicht genug: die Belichtung verändert auch die Signalvarianz. Gegeben ein weißes Signal $X$, wird seine lokale Varianz $V_1$ über ein Abtastfenster $\mathcal{N}$ ausgedrückt 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} $$

Wenn wir statt $X$ dasselbe Bild um einen Faktor $l$ überbelichtet hätten, dann würde die Varianz von $lX$ zu $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} $$

Die Varianz des Signals steigt also mit dem Quadrat des Belichtungsfaktors. Das ist hier aus zwei Gründen für uns von Bedeutung, die sich bereits jetzt so zusammenfassen lassen: alles, was wir hier tun, ist, die skalenweise Varianz des Bildes zu verändern.

Erstens ist unser B-Spline-Unschärfeschritt ein lokaler gewichteter Mittelwert, und $H_s$, ausgewertet am Pixel mit den Koordinaten $(x, y)$, kann tatsächlich geschrieben werden 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} $$

mit $h_{i, j}$ den Koeffizienten des 5-Tap-2D-B-Spline-Kerns und $X$ dem geglätteten Signal der vorherigen Skala (oder dem Ausgangsbild beim ersten Schritt). Wir können analog zeigen, dass $H_s$ linear vom Belichtungsfaktor abhängt. Die obige Gleichung zeigt, wie $H_s$ als Modulation um einen lokalen Mittelwert aufgefasst werden kann: der Betrag dieser Modulation ist nicht unabhängig vom Betrag des Signals. Das bedeutet, dass jede Glättung von $H_s$ (durchgeführt als diffusiver Prozess) je nachdem, ob das Bild über- oder unterbelichtet ist, ein anderes Gewicht und eine andere Auswirkung auf Details hat, obwohl der Inhalt derselbe ist.

Anders gesagt: erst glätten (oder umgekehrt schärfen) und dann unterbelichten, oder erst unterbelichten und dann glätten, wird nicht denselben Effekt auf Details haben, obwohl der endgültige Gesamtbetrag (Mittelwert) des Signals derselbe sein wird. Das ist nicht das, was wir von einem gutartigen Bildfilter erwarten: die Datenrepräsentation des Inhalts sollte nicht beeinflussen, wie wir den Inhalt selbst verarbeiten. In der diffusiven Situation ist das nicht so schädlich, aber in der Schärfungssituation werden Details in den Schatten im Vergleich zu Details in den Lichtern ohne jegliche Normalisierung wirklich überschärft.

Zweitens verändert die B-Spline-Unschärfe (oder ihre beste Gaußsche Näherung) auch die Varianz des Signals. Wenn wir das diskrete Signal $X$ als lokale Modulation um den globalen Mittelwert $\mu$ ausdrücken, erhalten wir $X_n = \mu + \epsilon_n$. Dann wird die auf $X$ angewandte B-Spline-Unschärfe zu:

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

Da die Koeffizienten des Kerns $h_{k}$ normalisiert sind und $\mu$ per Definition über das Fenster der Länge $k$ konstant ist, erhalten wir $\mu \sum_{k} h_{k} = \mu$, was bedeutet, dass die Unschärfe den Mittelwert nicht verändert. Die Varianz wird dann ausgedrückt 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} $$

Daraus können wir zeigen, dass für weiße, unkorrelierte Signale $\operatorname{Var}(G_0) = \operatorname{Var}(X) \sum_k h_k^2$ gilt, und allgemeiner:

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

wobei $\sum_k h_k^2 = (35 / 128)^2$ für den kardinalen 5-Tap-2D-B-Spline-Filter gilt.

Und auch hier haben wir ein Problem: dasselbe Objekt, abgetastet bei einer gewissen Auflösung oder bei der vierfachen dieser Auflösung, würde eine Wavelet-Zerlegungsstufe später bei derselben „Frequenz" erscheinen, wodurch seine Flächenvarianz um $(35 / 128)^2$ geringer wäre. Aber das ist rein ein Abtastartefakt. Die Varianz des Objekts selbst kann außerhalb des Bildes, in einer kontinuierlichen Betrachtung, als die lokalen Farbmodulationen um seine mittlere Oberflächenfarbe konzeptualisiert werden. Auch wenn diese idealisierte Varianz von keinem Abbildungsapparat wiederhergestellt werden kann, sollte die Art, wie wir das Signal behandeln, zumindest in der Varianz stabil sein, damit wir postulieren können, dass die Bildvarianz die idealisierte Varianz des Objekts bis auf einen konstanten Skalierungsfaktor repräsentiert.

Das mag wie eine philosophische Sorge erscheinen, bis wir auf ein praktisches Problem jeder Bildbearbeitungssoftware stoßen: was passiert, wenn wir den Effekt hinein-/herausgezoomt in der Vorschau betrachten? Wie skalieren wir den Effekt, damit die herunterskalierte Vorschau dem Ergebnis in voller Auflösung noch treu bleibt?

All diese Abtastungsdiskrepanzen müssen also normalisiert werden, um einen Bildfilter zu erreichen, der versucht, Inhalt unabhängig von seiner Datenrepräsentation zu manipulieren.

Definition einer Regularisierungsmetrik

Wir haben oben gesehen, wie die Signalvarianz eine relevante Metrik für das ist, was wir hier tun: wir können sie durch die Unschärfeschritte hindurch verfolgen, sie mit dem Signalbetrag verknüpfen, und sie repräsentiert die Signalmodulation um den Mittelwert.

Leider haben wir keinen Zugang zu einer Varianzmetrik, sobald wir in das Schema der Wavelet-Zerlegung eintreten. Wir haben jedoch oben gesehen, dass $H_s$ konzeptionell dem Term $(\bar{X} - X_i)$ der Varianz recht nahekam:

• statt eines arithmetischen Mittels verwenden wir einen gewichteten Mittelwert unter Verwendung von B-Spline-Koeffizienten,
• statt eines globalen Mittelwerts verwenden wir einen lokalen,
• die radiale Natur des B-Splines macht ihn rotationsinvarianter als jeder quadratische patchweise Mittelwert.

Wir werden daher die Energie des Bandes $H_s$ verwenden, ausgewertet an denselben Pixelkoordinaten wie das Laplace-Stencil, definiert als:

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

Für ein langsam variierendes, weißes, unkorreliertes Signal kommt $\overline{Q_s} = Q_s / |\mathcal{N}_{3\times 3}|$ der skalen- und patchweisen Varianz nahe.

Die Regularisierung ist für das Schärfungsproblem gedacht, das schlecht gestellt ist: in dieser Situation erhöhen wir die Energie jeder Schicht $H_s$ und benötigen einen Parameter, um sie irgendwann in Schach zu halten. Das ist ein übliches Vorgehen bei inversen Problemen wie Entrauschen und Entschärfen (deblurring), für die die Totalvariation (Total Variation) schon geraume Zeit als Regularisierungsschema verwendet wird.

Das Regularisierungsmodell, das wir verwenden werden, ist:

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

mit den Benutzerparametern:

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

Wir haben oben gezeigt, wie die Signalvarianz mit dem Quadrat der Belichtungsskalierung variiert und wie $H_s$ linear mit der Belichtungsskalierung variiert. $G_s$ trägt durch seine Eigenschaft, ein lokaler gewichteter Mittelwert zu sein, dieselbe lineare Abhängigkeit.

Das Verhältnis $H_s / G_s$ ist also belichtungsinvariant. Durch Identifikation gilt $L_s(q) = G_s(q)$ in der Regularisierungsgleichung, sodass wir die belichtungsinvariante Bandenergie verwenden:

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

und ihren lokalen Mittelwert:

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

Normalisierung von Skala und räumlicher Abdeckung

Das 3×3-Laplace-Stencil dehnt sich für jede Skala $s$ um eine Schrittweite von $2^s$ aus, genau wie der B-Spline-Kern: das ist die Grundlage des „à-trous"-Schemas. Der von diesem Kern abgedeckte physische Raum wächst mit den Wavelet-Skalen.

Wir können zeigen, dass ein 2D-Gaußfilter mit Parameter $\sigma$ die ursprüngliche Signalvarianz wie folgt verändert:

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

Da $4 \pi \sigma^2$ die effektive Scheibenfläche ist, die vom Gaußfilter abgedeckt wird, passt das zu einer sehr einfachen Intuition: der Gaußfilter breitet die ursprüngliche Modulation (ausgedrückt als Varianz) über eine größere Fläche aus. Das ist Diffusion in aller Kürze. Daraus leiten wir ab:

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

Zwischen 2 Gaußschen Unschärfeschritten variiert der Varianzparameter des äquivalenten Gaußfilters um:

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

Und von einer Skala zur anderen wächst dieser Radius um $\Delta \sigma = \sigma_B 2^s \approx 2^s$, unter Berücksichtigung von $\sigma_B \approx 1.05\dots$.

Die Auswertung des Laplace dehnt sich also räumlich mit derselben Rate aus (dieselbe Schrittweite bei jeder Skala) wie die B-Spline- oder Gaußfilter: der Laplace folgt implizit der Ausbreitung der Varianz über die Skalen hinweg, und hier ist keine zusätzliche Skalennormalisierung nötig.

Die anfängliche Implementierung von diffuse or sharpen hatte einen $\sigma^2_s$-Boost, der auf den Regularisierungsparameter $\lambda$ angewandt wurde: in der Praxis verhindert das, dass grobe Skalen überhaupt sichtbare Auswirkungen haben, und beschränkt sie darauf, die Grundlage für hochfrequente Arbeit zu sein.

Ein weiterer Versuch wurde unternommen, $\lambda$ zu normalisieren, unter Berücksichtigung der Tatsache, dass $Q_s$ eine mit steigendem $s$ abnehmende Bandenergie hat, sodass eine skaleninvariante Energiemetrik zu Folgendem führte:

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

Selbst kombiniert mit dem $\sigma^2_s$-Faktor oben ergab das viel zu viel Gewicht auf den groben Skalen, was es schwierig machte, bei hohen Frequenzen zu schärfen, während die niedrigen bereits stark überschossen (was einen verschmierten Eindruck ergab). Beide dieser Versuche wurden aufgegeben.

Ergebnis

After

Before

Das obige Bild wird den meisten übertrieben erscheinen, aber das ist nicht der Punkt: ich habe meinen Anteil an kaputten Schärfungsoperatoren gehabt, die gut genug funktionierten, solange man die Stärke nicht über 2 % trieb. Solche Bildfilter sind völlig uninteressant. Man lernt viel mehr über Algorithmen, indem man betrachtet, wie sie versagen, als indem man betrachtet, wie sie in ihrem optimalen Bereich Erfolg haben.

Dieses Bild ist die ultimative Falle für jeden Schärfungsalgorithmus:

• der Vordergrund ist näher und weniger dunstig als der Hintergrund, sodass er geradezu darum bettelt, überschärft zu werden,
• der nahe Vordergrund ist zudem viel dunkler, sodass es dort ebenfalls leicht ist, die Schärfe bis in den Himmel hochzutreiben,
• wir haben einen sehr kontrastreichen Bergkamm, der geradezu darum bettelt, Halos um Kanten zu erzeugen,
• die Sonnenscheibe in einem bewölkten Himmel wird von den meisten Algorithmen mit dunklen Kanten geschärft,
• der erforderliche Grad an Entdunstung würde das Rauschen explodieren lassen (dies wurde 2017 mit einer Nikon D90 bei 200 ISO aufgenommen, wir sind weit von aktuellen Sensoren entfernt).

Hier führen wir also ein gemeinsames Entrauschen und Entschärfen (denoising & deblurring) bei großen Radien durch. Das Entschärfen nutzt isotrope Gegendiffusion (die den atmosphärischen Dunst am besten beschreibt), auf Hochfrequenz und Niederfrequenz gleichermaßen. Das Entrauschen nutzt Isophoten-Diffusion auf der Hochfrequenz, dem in der Hochfrequenz abgetasteten Gradienten folgend. All das wurde ohne Maskierung in einer einzigen diffuse or sharpen-Instanz durchgeführt.

Insgesamt sehen wir kein Kantenüberschießen und keine Halos. Der dunkle Vordergrund wurde größtenteils ignoriert, wie es sein soll, aber wir haben die Details unten im Tal zurückbekommen. Wir haben keine Farbverschiebung und keine chromatischen Aberrationen. Ich werde nicht behaupten, dass dieses Bild artefaktfrei ist, denn bei genauem Hinsehen bekommen wir Farbstreifen und neue Details, die entweder tatsächliche Rekonstruktion oder schlichte Halluzinationen des diffusiven Modells sein könnten. Aber es bleibt festzuhalten, dass diese Artefakte, falls vorhanden, organisch genug aussehen, um unbemerkt zu bleiben, sofern wir das Originalbild nicht haben.

Dies bestätigt die Relevanz des diffusiven Multiskalen-Schemas zusammen mit der Regularisierungsstrategie.

After

Before

After

Before

Dieses wurde mit einer Sony Alpha 7R Mk 2, aber einem Konika Hexanon 40 mm f/1.8 aufgenommen. Das ist ein Pancake-Objektiv aus den 1970er-Jahren, und so ist das Bild 45 Mpx wunderschöne Unschärfe. Porträts sind weniger nachsichtig als Landschaften, weil die Haut gesund bleiben muss, wenn man die Schärfe zurückbringt. Hier haben wir zudem all die Lensflares und die Sonnenstrahlen, die schnell entarten können.

Das Detail ist hier besonders beeindruckend, weil wir die Augenbrauen und Wimpern zurückbekommen, ohne das hinterleuchtete Haar zu überschärfen. Insgesamt gibt es kein Ringing, Fringing oder Haloing. Leider hat sich die Hautqualität verschlechtert, und weiter kommen wir nicht, ohne die Haut zu maskieren, um sie auszuschließen. Von da an geht es nur noch darum, mit einer Tonwertkurve selektiv den globalen Kontrast zurückzubringen:

After

Before

Was spezifisch für Ansel ist

Die wissenschaftlichen Referenzen erklären die Bausteine:

• anisotrope Diffusion und Wärmeübertragungs-Inpainting,²
• Skalenraum und diskrete Diffusion,³⁴
• isotrope 9-Punkt-Laplace-Operatoren,⁵⁶
• B-Spline-Gaußnäherung.⁷

Was spezifisch für Ansel ist, ist die Art, wie sie zusammengesetzt werden:

• eine à-trous-B-Spline-Pyramide in voller Auflösung,
• vier unabhängig parametrisierte Diffusionsoperatoren, aufgeteilt zwischen niedrigen und hohen Frequenzen,
• ein HF-Bandenergie-Regularisierer, zusätzlich normalisiert durch die lokale LF-Energie,
• zoom-bewusste Skalenauswahl in der Dunkelkammer-Vorschau,
• identische CPU- und OpenCL-Mathematik.¹

Das Modul sollte also als ingenieurmäßige Synthese mehrerer numerischer Ideen verstanden werden, nicht als wörtliche Umsetzung einer einzelnen Veröffentlichung.

Perspektiven

Das Multiskalen-Wavelet-Schema der anisotropen Diffusions-PDE mit Regularisierung liefert brauchbare fotografische Ergebnisse, über den bloßen Machbarkeitsnachweis hinaus. Es ermöglicht, gemeinsame Schärfung und Entrauschung zusammen mit regulärer orientierter Diffusion zu nutzen. Es erlaubt auch, per Software alte Objektive wiederzubeleben, die für die hochauflösende Digitalfotografie als ungeeignet galten. Das Modul selbst bietet einen generischen PDE-Spielplatz, der auf viele verschiedene Weisen genutzt werden kann.

Das Problem ist jedoch, dass die Natur der Einstellungen in (mindestens) Mathematik auf Bachelor-Niveau verwurzelt und für die meisten Fotografen kryptisch ist. Die Wirkung der Parameter zu erklären ist schwierig, ohne einzutauchen in das, was sie mathematisch bedeuten. Der Versuch, sie nach ihrer Funktion statt nach ihrer Natur umzubenennen, ist zum Scheitern verurteilt, weil ihre Funktion davon abhängt, wie sie miteinander kombiniert werden, und ob sie im positiven oder negativen Wertebereich verwendet werden.

Der Diffusionsaufbau ist ziemlich unkompliziert und erfordert keine Regularisierung. Er benötigt nicht alle 4 Ordnungen zusammen, sondern nur die Einstellungen erster und dritter Ordnung genügen. In der isotropen Situation ist er vollständig äquivalent zu einer Gaußschen Unschärfe, die (weil nicht-iterativ) weniger aufwendig zu berechnen sein wird.

Der Schärfungsaufbau ist zusammen mit dem gemeinsamen Entrauschen komplizierter. Der einzige Weg, ihn benutzerfreundlicher zu machen, besteht darin, ihn als Machine-Learning-Algorithmus zu trainieren:

• Paare von sauberen und unscharfen/verrauschten/dunstigen Bildern derselben Szene aufnehmen (Bewegungsunschärfe, Objektiv-Defokusunschärfe, weiche Objektive),
• einen Brute-Force-Parameter-Sweep des Moduls diffuse or sharpen durchführen, das die Rekonstruktion der schmutzigen Bilder vornimmt, und die $L_2$-Norm des Fehlers zwischen den sauberen Referenzbildern und den versuchten Rekonstruktionen aufzeichnen,
• die Bilder manuell in binäre Kategorien (unscharf, verrauscht, dunstig) oder nach Intensität klassifizieren (Rauschen kann als PSNR, RMS, … gemessen werden; der Radius wird wahrscheinlich auch dabei sein müssen),
• das Machine-Learning-Problem wird zu: was sind im 14D-Raum der Eingabeparameter von diffuse or sharpen die 4 Hauptrichtungen (verbunden mit Entrauschen, Entschärfen, Entdunsten, Radius), die den Fehler $E = ||\text{clean} - \text{reconstructed}||_2$ minimieren? Wir suchen nach den Eigenvektoren dieser Hauptrichtungen.
• das durch gewichtete PLS (Partial Least-Squares) lösen: $Y = X’ B + c$ für $n$ Paare sauberer/rekonstruierter Bilder, wobei:
  • $Y$ der $4 × n$-Vektor der Kategorien für jede Stichprobe ist,
  • $X’$ der $14 × n$-standardisierte Vektor (komponentenweise standardisiert: $X_j’ = \frac{X_r - \mu}{\sigma}$) der Parameter von diffuse or sharpen ist,
  • $B$ die $14 × 4$-Matrix ist, die unsere 14 kryptischen D or S-Parameter auf 4 benutzerfreundliche Parameter abbildet (das ist hier die Unbekannte),
  • $c$ das Residuum ist (Skalar oder Vektor, je nachdem, was passt),
  • die $L_2$-Norm des Fehlers als PLS-Gewichtung verwendet wird (wahrscheinlich in eine Exponentialfunktion eingespeist),
• sobald die Matrix $B$ bekannt ist (Abbildung 14D -> 4D), sie invertieren, um das Modell 4D -> 14D zu erhalten,
• einen alternativen GUI-Modus hinzufügen, der die 4 benutzerfreundlichen Parameter freilegt, sowie eine GUI-<->-Parameter-Schicht, die diese 4 in die 14 Eingabeargumente von D or S umwandelt (das heißt: das Matrixprodukt schreiben).

Jeder andere Versuch, diffuse or sharpen zu „vereinfachen", wird nur eine alberne Neubeschriftungsarbeit sein, die die tatsächliche Bedeutung der Parameter verschleiert und jeden mit dem entsprechenden mathematischen Hintergrund daran hindert, sie zu verstehen. Es wäre eine Schande, dafür zu sorgen, dass die wenigen Leute, die sie verstehen können, tatsächlich davon abgehalten würden, es überhaupt zu versuchen. Im Moment ist D or S schwer zu verstehen, aber immerhin kann es erklärt werden. Das Neubeschriften von Bedienelementen wird es nicht leichter verständlich machen, sondern nur eine zusätzliche Schicht semantischer Übersetzung zwischen Mathematik und GUI hinzufügen, die höchstwahrscheinlich ohnehin ungenau und irreführend ist, was es nur kognitiv anspruchsvoller machen wird, es zu erklären und zu erfassen.