Matematyka modułu diffuse or sharpen

Ten artykuł dokumentuje matematykę aktualnie zaimplementowaną w module diffuse or sharpen Ansela, w postaci, w jakiej występuje w plikach src/iop/diffuse.c, src/common/bspline.h, data/kernels/diffuse.cl oraz data/kernels/bspline.cl. Nie jest to poradnik użytkowania. To rekonstrukcja modelu numerycznego na podstawie kodu źródłowego, w której twierdzenia naukowe zostały prześledzone aż do odwołań cytowanych w komentarzach kodu.¹

Streszczenie

Moduł działa na wieloskalowym rozkładzie à-trous obrazu w pełnej rozdzielczości, zbudowanym z rozmyć kardynalnymi B-sklejanymi. Na każdej skali Ansel oblicza cztery anizotropowe operatory drugiego rzędu, osobno dla pasm niskoczęstotliwościowych i wysokoczęstotliwościowych, reguluje ich działanie za pomocą znormalizowanej względem skali lokalnej energii pasma wysokoczęstotliwościowego, opcjonalnie ogranicza aktualizację do binarnej maski uzupełniania (inpainting) i rekonstruuje obraz od skali zgrubnej do drobnej.² Dyskretne operatory przestrzenne to szablony różnic centralnych 3x3 próbkowane na rzadkiej podkracie à-trous bieżącej skali, więc większe fizyczne wsparcie wynika z kroku $2^s$ drabiny falkowej, a nie z większego jądra PDE.³¹

Motywacja

Celem modułu diffuse było zasymulowanie akwarelowego pigmentu wypływającego z obrazu na jego brzegi, o tak :

Okazało się, że sama zmiana znaku aktualizacji w równaniu pochodnych cząstkowych dyfuzji może w rzeczywistości wyostrzyć obraz. Przez lata pracowałem nad ślepą dekonwolucją, zanim zająłem się modelem dyfuzji, i nigdy nie udało mi się sprawić, by nawet metoda state-of-the-art działała poza ładnymi, plakatowymi przypadkami (niskie ISO, niskie nasycenie, jednorodne rozmycie). Ale to oczywiście wymagało parametru regularyzacji, aby uniknąć rozbieżności.

Tak właśnie narodził się moduł diffuse or sharpen, jako ogólny szkielet do symulacji wszelkiego rodzaju zjawisk dyfuzyjnych lub kontr-dyfuzyjnych.

Model ciągły

W swej istocie moduł implementuje rodzinę anizotropowych równań dyfuzji postaci

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

gdzie $u$ to obraz, a $\mathbf{A}$ to symetryczny dodatni tensor dyfuzji, który steruje wygładzaniem albo wzdłuż izofot, albo wzdłuż gradientu, albo izotropowo.²⁴

Pierwotną inspiracją w kodzie jest model uzupełniania (inpainting) Qina i in., którzy sprzęgają rekonstrukcję struktury i tekstury obrazu za pomocą anizotropowego transferu ciepła.² Ansel zachowuje tego ducha, ale stosuje PDE w dziedzinie falkowej, osobno na pasmach niskoczęstotliwościowych i wysokoczęstotliwościowych, oraz udostępnia cztery sterowane przez użytkownika współczynniki transportu.

Dyskretne gradienty i anizotropia

Dla każdego piksela i każdego kanału moduł wyodrębnia sąsiedztwo 3x3 i oblicza scentrowane pierwsze pochodne

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

Jest to standardowy gradient różnic centralnych używany w dyskretnej przestrzeni skali i metodach dyfuzji.³

Lokalna orientacja gradientu wynosi wówczas :

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

z typowym rozwiązaniem awaryjnym $(1,0)$, gdy moduł gradientu znika.¹ Zauważ, że unikniemy kosztownego obliczania kąta $\theta$ jako funkcji $\arctan2$ składowych $(x, y)$, ponieważ nigdy nie będzie on używany bezpośrednio w dalszej części.

Siła anizotropii jest przeliczana z parametru użytkownika $a$ na dodatni współczynnik

$$ \alpha = a^2, $$

podczas gdy znak $a$ wybiera tryb:

• $a = 0$: dyfuzja izotropowa,
• $a > 0$: dyfuzja zgodna z kierunkiem izofoty,
• $a < 0$: dyfuzja zgodna z kierunkiem gradientu.¹

Człon tłumiący wynosi wówczas

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

co odpowiada anizotropowej dyfuzji Qina i in. dla uzupełniania (inpainting).²

Dla dyfuzji zgodnej z izofotą tensor zapisany w lokalnej bazie obrazu wynosi :

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

Dla dyfuzji zgodnej z gradientem Ansel używa postaci odwróconej :

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

Szablony dyfuzji 3x3

W zależności od parametrów użytkownika, dla każdego piksela obliczymy szablon $\mathbf{K(A)}$ przy użyciu jednego z powyższych tensorów $\mathbf{A}$.

Mając symetryczny tensor $\mathbf{A}$ :

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

budujemy dyskretne, obrócone, anizotropowe jądro laplasjanu $\nabla \cdot (\mathbf{A}\nabla u)$ jako szablon 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}, $$

z dokładnością do konwencji znaku wspomnianej w kodzie źródłowym dla członów pozadiagonalnych.²⁴¹

W przypadku izotropowym rozwija się on do stałej :

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

co jest klasycznym izotropowym laplasjanem w stylu Oono-Puri. Jego zaletą jest zachowanie rotacyjne: w porównaniu z 5-punktowym laplasjanem błąd kątowy jest zredukowany, co ma znaczenie, gdy dyfuzja nie może faworyzować osi obrazu.⁵⁶ Ten błąd kątowy oceniono względem podobnych jąder w Rotation-invariant Laplacian for 2D grids .

Aby zrozumieć wpływ kierunku na dyfuzję, zacznijmy od zaszumionego dysku i poddajmy go dyfuzji za pomocą Ansela :

Warto tu wspomnieć o kilku rzeczach :

• Tylko rozmycie izotropowe daje achromatyczny dysk, pozostałe zachowują pewien wielkoskalowy szum chromatyczny, ponieważ szum chromatyczny tworzy lokalne zmiany gradientu.
• Dyfuzja anizotropowa jest doskonała w pobliżu pozycji północ/wschód/południe/zachód okręgu, to znaczy tam, gdzie kąt gradientu jest idealnie zgodny z siatką pikseli (0° lub 90°). Pomiędzy nimi widzimy rozbieżności w sposobie traktowania kierunku, wynikające z ograniczeń numerycznych obróconego kwadratowego szablonu 3×3.

Piramida B-sklejana à-trous

Analiza wieloskalowa jest zbudowana z separowalnego 5-odczepowego filtra kardynalnego B-sklejanego

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

Filtr ten jest zwartym przybliżeniem gaussowskim, dlatego pojawia się naturalnie w metodach przestrzeni skali opartych na sklejanych.⁷ Wykorzystujemy tu ponownie szkielet wprowadzony przez Johannesa Hanikę dla modułu contrast equalizer.⁸

Obecna implementacja wykorzystuje historyczną nieredukowaną drabinę à-trous:

1. Pierwszy dolnoprzepustowy filtr w pełnej rozdzielczości to

   $$ G_0 = h_0 * u, $$

2. Najdrobniejsze pasmo szczegółów to

   $$ H_0 = u - G_0. $$

3. Grubsze poziomy zachowują tę samą rozdzielczość obrazu i jedynie zwiększają krok rozmycia o $2^s$:

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

   co oznacza, że przy $s=1$ odczepy $h_{2}$ są obliczane co drugi piksel, przy $s=2$ odczepy $h_4$ są obliczane co 4 piksele i tak dalej…

4. Pasmo przechowywane na skali $s$ to różnica między dwoma kolejnymi filtrami dolnoprzepustowymi:

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

Tak więc każde pasmo istnieje na oryginalnej siatce obrazu (bez decymacji), a PDE rozwiązujemy od drobnej do zgrubnej skali równolegle, akumulując rozwiązanie skala po skali do bufora wyjściowego, aż dodamy końcową resztę.¹ Ten schemat unika nakładania korekcji drobnych skal na skale zgrubne, jak to jest w zwyczaju przy piramidach gaussowskich lub solverach wielosiatkowych, i wykazał empirycznie lepszą stabilność w odwrotnym problemie rekonstrukcji ostrości. Jest to równoważne osobnemu zarządzaniu energią każdego pasma, w sposób przypominający korektory (equalizery) audio HiFi, lecz w ramach dwuwymiarowego szkieletu przestrzennego.

Skala jądra, skala pasma i obwiednia skali w GUI

Aby zbadać właściwości schematu falkowego, wykorzystamy właściwość kardynalnej B-sklejanej polegającą na tym, że jest ona przybliżeniem filtra gaussowskiego. Parametr gaussowski $\sigma$ kontroluje sposób składania jąder rozmycia:

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

więc wariancje $\sigma^2$ gaussowskich jąder sumują się przy splocie.³⁷ Natomiast wariancja przefiltrowanego sygnału $X$ :

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

zależy od widma $X$ i zasadniczo nie jest równa $\sigma^2$. Aby ujednoznacznić obie wariancje, w dalszej części tego artykułu $\sigma$ odnosi się wyłącznie do parametru jądra gaussowskiego, nigdy do pierwiastka kwadratowego wariancji sygnału.

Aby uniknąć niejednoznaczności, pomocne jest również oddzielenie poziomów dolnoprzepustowych od pasm poziomów:

• $G_0 = h_1 * u$ to pierwszy dolnoprzepustowy filtr w pełnej rozdzielczości nałożony na obraz wejściowy $u$,
• $G_s = h_{2^s} * G_{s-1}$ dla $s \ge 1$,
• $H_0 = u - G_0$,
• $H_s = G_{s-1} - G_s$ dla $s \ge 1$.

Tak więc:

• indeks dolnoprzepustowy $s$ odnosi się do poziomu rozmycia $G_s$,
• indeks pasma $s$ odnosi się do pasma szczegółów $H_s$ leżącego między $G_{s-1}$ a $G_s$,
• każdy poziom pozostaje próbkowany na oryginalnej siatce obrazu.¹

Filtr B-sklejany najlepiej przybliża równoważne jądro gaussowskie o parametrze :

$$ \sigma_B \approx 1.055365. $$

Ponieważ wariancje jąder gaussowskich sumują się przy splocie, efektywny promień rozmycia $G_s$ wynika bezpośrednio z drabiny analizy à-trous. Dlatego skumulowana sekwencja dolnoprzepustowa wynosi

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

gdzie

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

oraz, ogólnie,

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

Równoważnie,

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

Kod używa $\sigma_{G_s}$, aby zdecydować, ile skal jest potrzebnych do dopasowania do żądanego przez użytkownika parametru radius oraz aby zważyć każde pasmo wokół wybranego przez użytkownika centralnego promienia:

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

gdzie $z$ to poziom powiększenia w ciemni (w porównaniu z rawem w pełnej rozdzielczości), $r_c$ to radius_center w GUI, a $r_w$ to radius.¹ Ustawiając te wartości, użytkownicy definiują zanik pewnego rodzaju dyskretyzowanego filtra „szerokości pasma" wyśrodkowanego na dowolnej częstotliwości. W tym sensie $r_w$ to szerokość gaussowskiej obwiedni w przestrzeni skali. Nie jest to promień szablonu PDE ani promień pojedynczego rozmycia. To rozpiętość profilu wzmocnienia stosowanego do dyskretnych pasm falkowych.

Schemat ważenia jest stabilny niezależnie od poziomu powiększenia : promienie są brane w przestrzeni obrazu raw (pełnej rozdzielczości). Przy podglądzie obrazów pomniejszonych w ciemni najwyższe częstotliwości są obcinane przez pomniejszanie, więc stosujemy rozkład falkowy zaczynając od najwyższej dostępnej częstotliwości i przeważamy je zgodnie z ich równoważnym promieniem w pełnej rozdzielczości. Pozwala to na dość dokładny pomniejszony podgląd, choć może on ukrywać zaszumione artefakty pojawiające się na najwyższych poziomach.

Pasma, których równoważny promień rozmycia $\sigma_{G_s}$ leży blisko $r_c / z$, otrzymują największe wzmocnienie, podczas gdy pasma bardziej oddalone są stopniowo tłumione. Małe $r_w$ daje wąski wybór skal; duże $r_w$ daje szerszą i bardziej płaską odpowiedź w sąsiednich pasmach.

Zatem parę (radius_center, radius) należy odczytywać jako środek i szerokość w przestrzeni skal. Moduł nie definiuje bezpośrednio funkcji Gaussa w częstotliwości Fouriera; zamiast tego definiuje gaussowską obwiednię nad dostępnymi pasmami wieloskalowymi.¹

Ogólna wieloskalowa aktualizacja PDE

Niech:

• $H_s$ będzie zapamiętanym pasmem szczegółów przy indeksie pasma $s$,
• $G_s$ będzie bieżącą rekonstrukcją niskoczęstotliwościową używaną przy rozwiązywaniu pasma $H_s$ podczas syntezy.

Dla danego pasma $s$ kod buduje cztery odpowiedzi dyfuzyjne:

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

gdzie:

• $K_{1,s}$ i $K_{2,s}$ są anizotropowymi maskami dyfuzji 3x3 stosowanymi do bieżącej rekonstrukcji niskoczęstotliwościowej $G_s$,
• $K_{3,s}$ i $K_{4,s}$ są analogicznymi maskami stosowanymi do pasma szczegółów $H_s$,
• $a_1$ do $a_4$ to 4 zdefiniowane przez użytkownika współczynniki tłumienia anizotropii (patrz wyżej),
• $p_1$ do $p_4$ to 4 zdefiniowane przez użytkownika współczynniki aktualizacji (transportu) PDE (patrz niżej),
• splot jest wykonywany kanał po kanale na rzadkiej podsieci à-trous o kroku $2^s$.¹

Wszystkie jądra laplasjanu $K$ są stosowane z tą samą odległością kroku, co rozmycie B-sklejane użyte do utworzenia skali falkowej $s$, nad którą są stosowane, czyli $2^s$. Zważywszy, że jądro rozmycia ma rozmiar 5×5, a jądra laplasjanu 3×3, oznacza to, że laplasjan pokrywa jedną czwartą powierzchni rozmycia.

Wszystkie te $K$ są próbą uczynienia diffuse or sharpen ogólnym wieloskalowym frameworkiem PDE, ponieważ będą one krzyżować informacje pomiędzy $G_s$ a $H_s$ :

Gradienty obliczane na $G_s$ z większym prawdopodobieństwem wskazują na prawdziwe szczegóły obrazu i mogą być użyte w kontekście wyostrzania do ignorowania szumu. Gradienty obliczane na $H_s$ są bardziej wrażliwe na szum i mogą być użyte w kontekście dyfuzyjnym do zmiękczania szumu. Oddzielamy warstwę, na której wnioskujemy o strukturze obrazu (kierunek gradientu), od warstwy, na której obliczamy aktualizację PDE (laplasjan), a więc dyfuzję.

Z braku lepszego określenia, te odpowiedzi $K_i$ są powiązane z parametrami GUI zwanymi „rzędami" od pierwszego do czwartego :

• speed, czyli współczynnik aktualizacji PDE,
• anisotropy, czyli współczynnik tłumienia $a$ z tensora anizotropii powyżej.

Mamy zatem 4 parametry speed, po jednym dla każdego rzędu :

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

I podobnie, 4 parametry anisotropy : $$ (a_1, a_2, a_3, a_4) = (\texttt{first}, \texttt{second}, \texttt{third}, \texttt{fourth}) $$

Znaczenie tego wszystkiego można przełożyć na język potoczny w następujący sposób :

1. dyfundujemy strukturę w kierunku struktury ($D_1$, pierwszy rząd),
2. dyfundujemy strukturę w kierunku tekstury ($D_2$, drugi rząd),
3. dyfundujemy teksturę w kierunku struktury ($D_3$, trzeci rząd),
4. dyfundujemy teksturę w kierunku tekstury ($D_4$, czwarty rząd).

Dowolny współczynnik $p$ ustawiony na 0 znosi dyfuzję, dowolny współczynnik $a$ ustawiony na 0 znosi anizotropię.

Aktualizacja PDE na skali $s$ wynosi zatem :

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

gdzie :

• $\kappa$ jest współczynnikiem dyskretyzacji, czyli $\frac14$ dla skończonych różnic centrowanych,
• $\nu_s$ jest parametrem regularyzacji, który omówimy w następnej sekcji,
• $w_s$ jest wagą skali zdefiniowaną w poprzedniej sekcji.

A ostateczna resynteza to po prostu :

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

Zatem, jeśli podsumujemy cały algorytm, przyjmując $u$ jako obraz początkowy, $u’$ jako końcowy, $n$ jako końcową liczbę skal :

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

Kilka uwag:

• $n$ nie jest parametrem użytkownika, lecz jest wyznaczane względem docelowego końcowego $\sigma$ żądanego przez parametry promienia ustawione przez użytkownika. Jako efekt uboczny dostosowuje się to do poziomu powiększenia.
• Rozmycie Gaussa jest samo w sobie dwuwymiarowym izotropowym rozwiązaniem równania ciepła : jest już dyfuzją.
• Zwiększanie promienia rozmycia jest równoważne pozwoleniu, aby dyfuzja przebiegała przez dłuższy czas : sygnał rozprzestrzenia się dalej.
• Dla $s > 0$ pasmo $H_s$ w rzeczywistości staje się różnicą funkcji Gaussa. Przy pewnym współczynniku korekty skalowania różnica funkcji Gaussa jest przybliżeniem laplasjanu funkcji Gaussa, który sam jest estymacją laplasjanu przy współczynniku skalowania $\sigma$.
• Ponowne zastosowanie (być może anizotropowego) laplasjanu na $H_s$ jest równoważne pochodnej cząstkowej czwartego rzędu (bi-laplasjan).

Struktura współczynnika korekcyjnego

Wyrażenie sterujące transportem PDE to

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

W tej sekcji zdefiniujemy współczynnik $\nu_s$.

Definiowanie dobrze zachowujących się filtrów częstotliwości dla fotografii

Fotografie są cyfrowymi reprodukcjami obrazu utajonego, powstałymi przez aparaturę naświetlającą (przysłona, czułość ISO matrycy, czas otwarcia migawki) oraz aparaturę dyskretyzacji (czyli próbkowania przestrzennego) (matryca filtrów barwnych, siatka pikseli). Są to artefakty technologii użytej do uchwycenia obrazu utajonego i nie dotyczą właściwego obrazu.

Niestety, heurystyki przechwytywania dotyczące naświetlenia i próbkowania wpływają na sposób, w jaki przetwarzamy obraz cyfrowy. Pokazałem na dyskach w sekcji maski dyfuzji 3×3, jak krawędzie ukośne zachowują się inaczej niż krawędzie zgodne z siatką (pionowe/poziome), mimo że wybrałem jądra najbardziej niezmiennicze względem obrotu : obrót treści obrazu (względem siatki pikseli) zmieni sposób obliczania gradientów dyskretnych. Konkretną implikacją jest tutaj to, że : obrócenie obrazu przed lub po diffuse or sharpen nie da tego samego rezultatu.

Ale to nie koniec : naświetlenie zmienia także wariancję sygnału. Dla białego sygnału $X$ jego lokalna wariancja $V_1$ w oknie próbkowania $\mathcal{N}$ wyraża się :

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

Gdybyśmy zamiast przechwytywania $X$ prześwietlili ten sam obraz o czynnik $l$, to wariancja $lX$ stałaby się $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} $$

Zatem wariancja sygnału rośnie z kwadratem czynnika naświetlenia. Ma to dla nas tutaj znaczenie z dwóch powodów, które można już teraz podsumować tak : wszystko, co tutaj robimy, sprowadza się do zmiany wariancji obrazu w poszczególnych skalach.

Po pierwsze, nasz krok rozmycia B-sklejanego jest lokalną średnią ważoną, a $H_s$ obliczone w pikselu o współrzędnych $(x, y)$ można faktycznie zapisać :

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

gdzie $h_{i, j}$ to współczynniki 5-punktowego, dwuwymiarowego jądra B-sklejanego, a $X$ to rozmyty sygnał z poprzedniej skali (lub obraz początkowy dla pierwszego kroku). Podobnie możemy pokazać, że $H_s$ jest liniowo zależne od czynnika naświetlenia. Powyższe równanie pokazuje, jak $H_s$ można postrzegać jako modulację wokół lokalnej średniej : amplituda tej modulacji nie jest niezależna od amplitudy sygnału. Oznacza to, że każde wygładzenie $H_s$ (realizowane jako proces dyfuzyjny) będzie miało inną wagę i inny wpływ na szczegóły w zależności od tego, czy obraz jest prześwietlony, czy niedoświetlony, mimo że treść jest ta sama.

Inaczej mówiąc, wygładzenie (lub przeciwnie, wyostrzenie), a następnie niedoświetlenie, albo niedoświetlenie, a następnie wygładzenie nie będą miały tego samego wpływu na szczegóły, mimo że końcowa ogólna amplituda (średnia) sygnału będzie taka sama. Nie jest to to, czego oczekujemy od dobrze zachowującego się filtra obrazu : reprezentacja danych treści nie powinna wpływać na to, jak przetwarzamy samą treść. W kontekście dyfuzyjnym nie jest to aż tak szkodliwe, ale w kontekście wyostrzania szczegóły w cieniach zostają naprawdę nadmiernie wyostrzone w porównaniu ze szczegółami w światłach, bez żadnej normalizacji.

Po drugie, rozmycie B-sklejane (lub jego najlepsze przybliżenie gaussowskie) zmieni także wariancję sygnału. Jeśli wyrazimy dyskretny sygnał $X$ jako lokalną modulację wokół globalnej średniej $\mu$, otrzymamy $X_n = \mu + \epsilon_n$. Wtedy rozmycie B-sklejane zastosowane do $X$ staje się :

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

Ponieważ współczynniki jądra $h_{k}$ są znormalizowane, a $\mu$ jest z definicji stałe w oknie o długości $k$, otrzymujemy $\mu \sum_{k} h_{k} = \mu$, co oznacza, że rozmycie nie zmienia wartości średniej. Wariancja wyraża się wtedy :

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

Stąd możemy pokazać, że dla białych, nieskorelowanych sygnałów $\operatorname{Var}(G_0) = \operatorname{Var}(X) \sum_k h_k^2$, a bardziej ogólnie :

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

gdzie $\sum_k h_k^2 = (35 / 128)^2$ dla 5-punktowego kardynalnego dwuwymiarowego filtra B-sklejanego.

I tutaj znów mamy problem : ten sam obiekt próbkowany z pewną rozdzielczością lub z 4-krotnością tej rozdzielczości pojawiłby się przy tej samej „częstotliwości" jeden krok dekompozycji falkowej później, co sprawiłoby, że jego wariancja powierzchni byłaby $(35 / 128)^2$ niższa. Ale jest to czysto artefakt próbkowania. Wariancję samego obiektu można skonceptualizować, poza obrazem, w ujęciu ciągłym, jako lokalne modulacje koloru wokół jego średniego koloru powierzchni. Choć tej wyidealizowanej wariancji nie da się odtworzyć żadną aparaturą obrazującą, sposób, w jaki traktujemy sygnał, powinien być przynajmniej stabilny w wariancji, abyśmy mogli założyć, że wariancja obrazu reprezentuje wyidealizowaną wariancję obiektu z dokładnością do pewnego stałego współczynnika skalowania.

Może się to wydawać rozważaniami filozoficznymi, dopóki nie natrafimy na praktyczny problem każdego oprogramowania do obrazowania : co się dzieje, gdy oglądamy efekt w powiększeniu/pomniejszeniu ? Jak skalować efekt, aby pomniejszony podgląd nadal wiernie oddawał rezultat w pełnej rozdzielczości ?

Zatem wszystkie te rozbieżności próbkowania muszą zostać znormalizowane, aby uzyskać filtr obrazu, który stara się manipulować treścią niezależnie od jej reprezentacji danych.

Definiowanie metryki regularyzacji

Widzieliśmy powyżej, że wariancja sygnału jest istotną metryką dla tego, co tu robimy : możemy prześledzić ją przez kroki rozmycia, powiązać ją z amplitudą sygnału i reprezentuje ona modulację sygnału wokół wartości średniej.

Niestety, nie mamy dostępu do metryki wariancji, gdy wejdziemy w schemat dekompozycji falkowej. Widzieliśmy jednak powyżej, że $H_s$ było koncepcyjnie dość bliskie członowi $(\bar{X} - X_i)$ wariancji :

• zamiast średniej arytmetycznej używamy średniej ważonej z użyciem współczynników B-sklejanych,
• zamiast średniej globalnej używamy lokalnej,
• radialna natura B-sklejanego czyni go bardziej niezmienniczym względem obrotu niż jakakolwiek średnia po kwadratowym wycinku.

Użyjemy więc energii pasma $H_s$, obliczanej przy tych samych współrzędnych piksela co maska laplasjanu, zdefiniowanej jako :

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

Dla wolno zmieniającego się, białego, nieskorelowanego sygnału $\overline{Q_s} = Q_s / |\mathcal{N}_{3\times 3}|$ staje się bliskie wariancji w poszczególnych skalach i wycinkach.

Regularyzacja jest przeznaczona dla problemu wyostrzania, który jest źle zdefiniowany : w tym kontekście zwiększamy energię każdej warstwy $H_s$ i potrzebujemy parametru, aby w pewnym momencie ją okiełznać. Jest to powszechna procedura w problemach odwrotnych, takich jak odszumianie i odblurowywanie, dla których wariacja całkowita (Total Variation) była od dłuższego czasu stosowana jako schemat regularyzacji.

Model regularyzacji, którego użyjemy, to :

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

z parametrami użytkownika :

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

Pokazaliśmy powyżej, jak wariancja sygnału zmienia się z kwadratem skalowania naświetlenia oraz jak $H_s$ zmienia się liniowo ze skalowaniem naświetlenia. $G_s$ niesie tę samą zależność liniową dzięki swojej właściwości bycia lokalną średnią ważoną.

Zatem stosunek $H_s / G_s$ jest niezmienniczy względem naświetlenia. Przez identyfikację $L_s(q) = G_s(q)$ w równaniu regularyzacji, więc używamy niezmienniczej względem naświetlenia energii pasma :

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

oraz jej lokalna średnia :

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

Normalizacja skali i pokrycia przestrzennego

Maska laplasjanu 3×3 rozszerza się o krok $2^s$ dla każdej skali $s$, tak samo jak jądro B-sklejane : jest to podstawa schematu „à-trous". Przestrzeń fizyczna pokrywana przez to jądro rośnie wraz ze skalami falkowymi.

Możemy pokazać, że dwuwymiarowy filtr Gaussa o parametrze $\sigma$ zmienia wariancję sygnału pierwotnego następująco :

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

Zważywszy, że $4 \pi \sigma^2$ jest efektywną powierzchnią dysku pokrywaną przez filtr Gaussa, odpowiada to bardzo prostej intuicji : filtr Gaussa rozprowadza pierwotną modulację (wyrażoną jako wariancja) po większej powierzchni. To jest dyfuzja w pigułce. Wyprowadzamy z tego :

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

Pomiędzy 2 krokami rozmycia Gaussa parametr wariancji równoważnego filtra Gaussa zmienia się o :

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

A z jednej skali na drugą ten promień rośnie o $\Delta \sigma = \sigma_B 2^s \approx 2^s$, przypominając, że $\sigma_B \approx 1.05\dots$.

Zatem obliczanie laplasjanu rozszerza się przestrzennie w tym samym tempie (ten sam krok przy każdej skali) co filtry B-sklejane lub Gaussa : laplasjan niejawnie podąża za rozprzestrzenianiem się wariancji między skalami i nie jest tu potrzebna żadna dodatkowa normalizacja skali.

Początkowa implementacja diffuse or sharpen miała wzmocnienie $\sigma^2_s$ stosowane do parametru regularyzacji $\lambda$ : w praktyce uniemożliwia to grubym skalom wywieranie jakiegokolwiek widocznego wpływu i ogranicza je do bycia podłożem dla pracy wysokoczęstotliwościowej.

Podjęto kolejną próbę znormalizowania $\lambda$ z uwzględnieniem faktu, że $Q_s$ ma malejącą energię pasma wraz ze wzrostem $s$, więc niezmiennicza względem skali metryka energii doprowadziła do :

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

Nawet w połączeniu z powyższym czynnikiem $\sigma^2_s$ dawało to zdecydowanie zbyt dużą wagę na grubych skalach, utrudniając wyostrzanie przy wysokich częstotliwościach, podczas gdy niskie już mocno przestrzeliwały (dając rozmazany wygląd). Obie te próby zostały porzucone.

Rezultat

After

Before

The picture above will look overcooked to most but it’s beyond the point : I have had my share of broken sharpening operators that worked well enough as long as you didn’t push the strength above 2%. Such image filters have zero interest. You learn about algorithms much more by looking at how they fail than at how they succeed in their sweet spot.

This image is the ultimate trap for every sharpening algorithms :

• the foreground is closer and less hazy than the background, so it’s begging to be oversharpened,
• the close foreground is also much darker, so again, easy to push sharpness up the sky there,
• we have a very contrasted mountain ridge begging to produce halos around edges,
• the sun disc in a cloudy sky will be sharpened by most algorithms with dark edges,
• the amount of dehazing required would make the noise explode (this was taken in 2017 with a Nikon D90 at 200 ISO, we are far from current sensors).

So, here we perform a joint denoising & deblurring at large radii. The deblurring uses isotropic counter-diffusion (which best describes atmospheric hazing), on high-frequency and low frequency alike. The denoising uses isophote diffusion on high-frequency, following the gradient sampled in high-frequency. All that was done with no masking, in a single diffuse or sharpen instance.

Overall, we see no edge overshooting and no halos. The dark foreground was mostly ignored, as it should, but we got the details down in the valley back. We got no color shift or chromatic aberrations. I won’t claim this picture is artifact-free, because, when looking up close, we get color streaks and new details that might be actual reconstruction or plain hallucinations of the diffusive model. But the point remain that these artificats, if any, look organic enough to go unnoticed if we don’t have the original image.

This validates the relevance of the multi-scale diffusive scheme, along with the regularization strategy.

After

Before

After

Before

This one was taken with a Sony Alpha 7R Mk 2 but a Konika Hexanon 40 mm f/1.8. This is a pancake lens from the 1970’s, and so the picture is 45 Mpx of beautiful blur. Portraits are less forgiving than landscapes because skin needs to stay healthy when bringing back sharpness. Here we also have all the lensflare and the sun rays that can quickly degenerate.

The detail is especially impressive here because we get back the eyebrowes and eyelashes without oversharpening the backlit hair. Overall, there is no ringing, fringing or haloing. Unfortunately, the skin quality has degraded and it’s as far as we can go without masking skin to exclude it. From there, it’s only about bringing back global contrast selectively with a tone curve :

After

Before

What is specific to Ansel

The scientific references explain the building blocks:

• anisotropic diffusion and heat-transfer inpainting,²
• scale-space and discrete diffusion,³⁴
• isotropic 9-point Laplacians,⁵⁶
• B-spline Gaussian approximation.⁷

What is specific to Ansel is the way they are assembled:

• a full-resolution a-trous B-spline pyramid,
• four independently parameterized diffusion operators, split between low and high frequencies,
• an HF band-energy regularizer additionally normalized by the local LF energy,
• zoom-aware scale selection in darkroom preview,
• identical CPU and OpenCL mathematics.¹

So the module should be understood as an engineering synthesis of several numerical ideas, not as a literal implementation of a single paper.

Perspectives

The multiscale wavelets scheme of anisotropic diffusion PDE with regularization produces usable photographic results, beyond the mere proof of concept. It allows to leverage joint sharpening and denoising, along with regular oriented diffusion. It also allows to rejuvenate, through software, old lenses that were deemed unfit for high-definition digital photography. The module itself provides a generic PDE playground that can be used in many different ways.

The problem though is that the nature of the settings is grounded into (at least) bachelor-level mathematics and is cryptic to most photographers. Explaining the effect of the parameters is difficult without diving into what they mean mathematically. Trying to rename them by their function instead of their nature is doomed to fail because their function depends on how they are combined together, and whether they are used in the positive or negative value range.

The diffusion setup is pretty straightforward and doesn’t require regularization. It doesn’t need the 4 orders altogether, but only the first and third order settings are enough. In the isotropic setting, it is fully equivalent to a Gaussian blur, which will be less expensive to compute (because non-iterative).

The sharpening setup, along with joint denoising, is more complicated. The only way to make it more user-friendly is to train it as a machine-learning algorithm :

• shoot pairs of clean and blurry/noisy/hazy pictures of the same scene (motion blur, lens defocus blur, soft lenses),
• perform a brute-force parameters sweep of diffuse or sharpen module performing reconstruction of the dirty images, and record the $L_2$ norm of the error between the clean reference images and the attempted reconstructions,
• manually classify the images in binary categories (blurry, noisy, hazy) or by intensity (noise can be measured as PSNR, RMS, … ; radius will probably need to be there too),
• the machine-learning problem becomes : in the 14D space of diffuse or sharpen input parameters, what are the 4 main directions (associated with denoising, deblurring, dehazing, radius) that minimize the error $E = ||\text{clean} - \text{reconstructed}||_2$ ? We are looking for the eigenvectors of those main directions.
• solve that by weighted PLS (partial least-squares) : $Y = X’ B + c$ for $n$ pairs of clean/reconstructed images, where :
  • $Y$ is the $4 × n$ vector of categories for each sample,
  • $X’$ is the $14 × n$ standardized vector (standardized component-wise : $X_j’ = \frac{X_r - \mu}{\sigma}$) of diffuse or sharpen parameters,
  • $B$ is the $14 × 4$ matrix mapping our 14 cryptic D or S parameters to 4 user-friendly parameters (it’s the unknown here),
  • $c$ is the residual (scalar or vector depending what fits),
  • the $L_2$ norm of the error is used as the PLS weighting (probably injected into an exponential),
• once the matrix $B$ is known (mapping 14D -> 4D), invert it to get the 4D -> 14 D model,
• add an alternative GUI mode exposing the 4 user-friendly parameters, and a GUI <-> parameters layer converting those 4 to the 14 input arguments of D or S (meaning : write the matrix product).

Any other attempt at “simplifying” diffuse or sharpen will only be a silly relabelling job obfuscating the real meaning of the parameters and preventing anybody with the proper mathematical background from understanding it. It would be a shame to make sure that the only few people able to understand it would actually be deterred of even trying. Right now, D or S is difficult to understand, but at least it can be explained. Relabelling controls will not make it easier to understand, but will only add an extra layer of semantic translation between math and GUI, most likely inaccurate and misleading anyway, which will only make it more cognitively demanding to explain and to grasp.