Les mathématiques de diffusion ou netteté

Cet article documente les mathématiques actuellement implémentées dans le module diffuse or sharpen d’Ansel, telles qu’on les trouve dans src/iop/diffuse.c, src/common/bspline.h, data/kernels/diffuse.cl et data/kernels/bspline.cl. Il ne s’agit pas d’un guide d’utilisation. C’est une reconstruction du modèle numérique à partir du code source, avec les affirmations scientifiques rattachées aux références citées dans les commentaires du code.¹

Résumé

Le module repose sur une décomposition multi-échelle a-trous de l’image à pleine résolution, construite à partir de flous B-spline cardinaux. À chaque échelle, Ansel calcule séparément quatre opérateurs anisotropes du second ordre sur les bandes de basse et de haute fréquence, régularise leur action par une énergie locale de bande haute fréquence normalisée par l’échelle, restreint éventuellement la mise à jour à un masque binaire de retouche, puis reconstruit l’image du grossier vers le fin.² Les opérateurs spatiaux discrets sont des stencils 3x3 à différences centrées échantillonnés sur la sous-grille a-trous clairsemée de l’échelle courante ; le support physique plus large provient donc du pas $2^s$ de l’échellettes plutôt que d’un noyau d’EDP plus grand.³¹

Motivation

L’objectif initial du module de diffusion était de simuler un pigment de type aquarelle s’échappant d’une image vers ses bords, comme ceci :

Il s’est avéré que le simple fait de changer le signe de la mise à jour de l’équation aux dérivées partielles de diffusion pouvait en réalité accentuer la netteté de l’image. J’avais travaillé pendant des années sur la déconvolution aveugle lorsque j’ai commencé à travailler sur un modèle de diffusion, sans jamais réussir à faire fonctionner correctement, hors des beaux cas d’école, même la méthode state-of-the-art : faible ISO, faible saturation, flou uniforme. Mais cela nécessitait bien sûr un paramètre de régularisation pour éviter la divergence.

C’est ainsi qu’est né le module diffuse or sharpen, comme cadre générique de simulation de toutes sortes de phénomènes diffusifs ou contre-diffusifs.

Modèle continu

Au coeur du module, on trouve une famille d’équations de diffusion anisotrope de la forme

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

où $u$ est l’image et $\mathbf{A}$ un tenseur de diffusion symétrique positif qui oriente le lissage soit le long des isophotes, soit le long du gradient, soit de manière isotrope.²⁴

L’inspiration d’origine dans le code est le modèle de retouche de Qin et al., qui couple structure de l’image et restauration de texture au moyen d’un transfert thermique anisotrope.² Ansel en conserve l’esprit, mais applique l’EDP dans un domaine en ondelettes, séparément sur les bandes de basse et de haute fréquence, et expose quatre coefficients de transport contrôlés par l’utilisateur.

Gradients discrets et anisotropie

Pour chaque pixel et chaque canal, le module extrait un voisinage 3x3 et évalue les dérivées premières centrées

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

Il s’agit du gradient standard à différences centrées utilisé dans les méthodes discrètes d’espace d’échelle et de diffusion.³

L’orientation locale du gradient est alors :

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

avec le repli habituel vers $(1,0)$ lorsque la norme du gradient s’annule.¹ Notons que nous éviterons de calculer l’angle coûteux $\theta$ comme fonction $\arctan2$ des composantes $(x, y)$, puisqu’il ne sera jamais utilisé directement dans la suite.

L’intensité de l’anisotropie est convertie à partir du paramètre utilisateur $a$ en un coefficient positif

$$ \alpha = a^2, $$

tandis que le signe de $a$ sélectionne le mode :

• $a = 0$ : diffusion isotrope,
• $a > 0$ : diffusion alignée avec la direction des isophotes,
• $a < 0$ : diffusion alignée avec la direction du gradient.¹

Le terme d’amortissement est alors

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

ce qui correspond à la diffusion anisotrope de Qin et al. pour la retouche.²

Pour une diffusion alignée sur les isophotes, le tenseur écrit dans la base locale de l’image est :

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

Pour une diffusion alignée sur le gradient, Ansel utilise la forme inversée :

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

Les stencils de diffusion 3x3

Selon les paramètres de l’utilisateur, nous évaluerons pour chaque pixel le stencil $\mathbf{K(A)}$ en utilisant l’un des tenseurs $\mathbf{A}$ ci-dessus.

Étant donné le tenseur symétrique $\mathbf{A}$ :

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

nous construisons le noyau laplacien discret, tourné et anisotrope $\nabla \cdot (\mathbf{A}\nabla u)$ sous forme du stencil 3x3 :

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

à la convention de signe près mentionnée dans le code source pour les termes hors diagonale.²⁴¹

Dans le cas isotrope, il se réduit à une constante :

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

qui est le laplacien isotrope classique de type Oono-Puri. Son intérêt réside dans son comportement rotationnel : comparé au laplacien à 5 points, l’erreur angulaire est réduite, ce qui importe lorsque la diffusion ne doit pas privilégier les axes de l’image.⁵⁶ Cette erreur angulaire a été évaluée face à des noyaux similaires dans Rotation-invariant Laplacian for 2D grids .

Pour comprendre l’effet de la direction dans la diffusion, commençons par un disque bruité et diffusons-le avec Ansel :

Deux ou trois choses méritent d’être mentionnées ici :

• Seul le flou isotrope produit un disque achromatique, les autres conservent un bruit de chroma à grande échelle, car le bruit de chroma crée des variations locales de gradient.
• La diffusion anisotrope est parfaite près des positions nord/est/sud/ouest du cercle, c’est-à-dire lorsque l’angle du gradient est parfaitement aligné sur la grille de pixels (0° ou 90°). Entre les deux, on observe des écarts dans la manière dont la direction est traitée, dus aux limitations numériques d’un stencil carré 3×3 tourné.

La pyramide B-spline à-trous

L’analyse multi-échelle est construite à partir du filtre B-spline cardinal séparable à 5 coefficients

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

Ce filtre est une approximation gaussienne compacte, ce qui explique pourquoi il apparaît naturellement dans les méthodes d’espace-échelle fondées sur les splines.⁷ Nous réutilisons ici le cadre introduit par Johannes Hanika pour le module contrast equalizer.⁸

L’implémentation actuelle utilise l’échelle historique à-trous non décimée :

1. Le premier passe-bas à pleine résolution est

   $$ G_0 = h_0 * u, $$

2. La bande de détail la plus fine est

   $$ H_0 = u - G_0. $$

3. Les niveaux plus grossiers conservent la même résolution d’image et ne font qu’agrandir le pas du flou d’un facteur $2^s$ :

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

   ce qui signifie qu’à $s=1$, les coefficients de $h_{2}$ sont évalués un pixel sur deux, qu’à $s=2$, les coefficients de $h_4$ sont évalués tous les 4 pixels, et ainsi de suite…

4. La bande stockée à l’échelle $s$ est la différence entre deux passe-bas successifs :

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

Ainsi chaque bande vit sur la grille d’image d’origine (pas de décimation), et nous résolvons l’EDP du fin vers le grossier en parallèle, en accumulant la solution échelle par échelle dans le tampon de sortie jusqu’à ce que nous ajoutions le résidu final.¹ Ce schéma évite d’empiler la correction des échelles fines par-dessus les échelles grossières, comme c’est l’usage avec les pyramides gaussiennes ou les solveurs multigrilles, et a montré empiriquement une meilleure stabilité dans le problème inverse de reconstruction de la netteté. Cela équivaut à gérer séparément l’énergie de chaque bande, d’une manière qui rappelle les égaliseurs audio haute-fidélité, mais dans un cadre spatial 2D.

Échelle du noyau, échelle de la bande, et l’enveloppe d’échelle de l’interface

Pour étudier les propriétés du schéma d’ondelettes, nous utiliserons la propriété de la B-spline cardinale d’être une approximation de filtre gaussien. Le paramètre gaussien $\sigma$ contrôle la manière dont les noyaux de flou se composent :

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

de sorte que les variances $\sigma^2$ du noyau gaussien s’additionnent sous la convolution.³⁷ Par contraste, la variance d’un signal filtré $X$ :

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

dépend du spectre de $X$ et n’est généralement pas égale à $\sigma^2$. Pour lever l’ambiguïté entre les deux variances, dans le reste de cet article, $\sigma$ ne désigne que le paramètre du noyau gaussien, jamais la racine carrée de la variance du signal.

Pour éviter toute ambiguïté, il est également utile de séparer les niveaux passe-bas des bandes des niveaux :

• $G_0 = h_1 * u$ est le premier passe-bas à pleine résolution sur l’image d’entrée $u$,
• $G_s = h_{2^s} * G_{s-1}$ pour $s \ge 1$,
• $H_0 = u - G_0$,
• $H_s = G_{s-1} - G_s$ pour $s \ge 1$.

Ainsi :

• l’indice de passe-bas $s$ désigne le niveau de flou $G_s$,
• l’indice de bande $s$ désigne la bande de détail $H_s$ située entre $G_{s-1}$ et $G_s$,
• chaque niveau reste échantillonné sur la grille d’image d’origine.¹

Le filtre B-spline approxime au mieux le noyau gaussien équivalent de paramètre :

$$ \sigma_B \approx 1.055365. $$

Comme les variances du noyau gaussien s’additionnent sous la convolution, le rayon de flou effectif de $G_s$ découle directement de l’échelle d’analyse à-trous. Par conséquent, la séquence cumulée de passe-bas est

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

avec

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

et, de manière générale,

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

De manière équivalente,

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

Le code utilise $\sigma_{G_s}$ pour décider combien d’échelles sont nécessaires pour correspondre au paramètre radius demandé par l’utilisateur et pour pondérer chaque bande autour du rayon central sélectionné par l’utilisateur :

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

où $z$ est le niveau de zoom en chambre noire (comparé au raw à pleine résolution), $r_c$ est radius_center dans l’interface, et $r_w$ est radius.¹ En les réglant, les utilisateurs définissent la décroissance d’une sorte de filtre de « bande passante » discrétisé centré sur une fréquence arbitraire. En ce sens, $r_w$ est la largeur d’une enveloppe gaussienne dans l’espace-échelle. Ce n’est pas le rayon du stencil de l’EDP, ni le rayon d’un flou unique. C’est l’étalement du profil de gain appliqué aux bandes d’ondelettes discrètes.

Le schéma de pondération est stable quel que soit le niveau de zoom : les rayons sont pris dans l’espace de l’image raw (pleine résolution). Lors de la prévisualisation d’images sous-échantillonnées en chambre noire, les fréquences les plus hautes sont tronquées par le sous-échantillonnage, aussi n’appliquons-nous la décomposition en ondelettes qu’à partir de la fréquence la plus haute disponible, en les repondérant selon leur rayon équivalent à pleine résolution. Cela permet une prévisualisation sous-échantillonnée assez précise, même si elle peut masquer des artefacts bruités apparaissant aux niveaux les plus hauts.

Les bandes dont le rayon de flou équivalent $\sigma_{G_s}$ est proche de $r_c / z$ reçoivent le gain le plus élevé, tandis que les bandes plus éloignées sont progressivement atténuées. Un petit $r_w$ donne une sélection étroite d’échelles ; un grand $r_w$ donne une réponse plus large et plus plate sur les bandes voisines.

Ainsi le couple (radius_center, radius) doit se lire comme un centre et une largeur dans l’espace-échelle. Le module ne définit pas directement une gaussienne en fréquence de Fourier ; il définit plutôt une enveloppe gaussienne sur les bandes multi-échelles disponibles.¹

Mise à jour EDP multi-échelle générique

Soit :

• $H_s$ la bande de détail stockée à l’indice de bande $s$,
• $G_s$ la reconstruction basse fréquence courante utilisée lors de la résolution de la bande $H_s$ pendant la synthèse.

À une bande donnée $s$, le code construit quatre réponses de diffusion :

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

où :

• $K_{1,s}$ et $K_{2,s}$ sont des stencils de diffusion anisotrope 3x3 appliqués à la reconstruction basse fréquence courante $G_s$,
• $K_{3,s}$ et $K_{4,s}$ sont les stencils analogues appliqués à la bande de détail $H_s$,
• $a_1$ à $a_4$ sont 4 coefficients d’amortissement d’anisotropie définis par l’utilisateur (voir ci-dessus),
• $p_1$ à $p_4$ sont 4 coefficients de mise à jour EDP (transport) définis par l’utilisateur (voir ci-dessous),
• la convolution s’effectue canal par canal sur le sous-réseau à-trous creux de pas $2^s$.¹

Tous les noyaux laplaciens $K$ sont appliqués avec la même distance de pas que le flou B-spline ayant servi à produire l’échelle d’ondelette $s$ sur laquelle ils sont appliqués, c’est-à-dire $2^s$. Étant donné que le noyau de flou est 5×5 et que les noyaux laplaciens sont 3×3, cela signifie que le laplacien couvre un quart de la surface du flou.

Tous ces $K$ visent à faire de diffuse or sharpen un cadre EDP multi-échelle générique, car ils vont croiser l’information entre $G_s$ et $H_s$ :

Les gradients évalués sur $G_s$ sont plus susceptibles de pointer vers des détails légitimes de l’image, et peuvent être utilisés dans un contexte d’accentuation pour ignorer le bruit. Les gradients évalués sur $H_s$ sont plus sensibles au bruit et peuvent être utilisés dans un contexte de diffusion pour adoucir le bruit. Nous séparons la couche sur laquelle nous inférons la structure de l’image (direction du gradient) de la couche où nous calculons la mise à jour EDP (laplacien) et donc la diffusion.

Faute de meilleur terme, ces réponses $K_i$ sont liées à des paramètres de l’interface appelés « ordres » du premier au quatrième :

• speed, qui est le coefficient de mise à jour EDP,
• anisotropy, qui est le coefficient d’amortissement $a$ du tenseur d’anisotropie ci-dessus.

Nous avons donc 4 paramètres speed, un pour chaque ordre :

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

Et de même, 4 paramètres anisotropy : $$ (a_1, a_2, a_3, a_4) = (\texttt{first}, \texttt{second}, \texttt{third}, \texttt{fourth}) $$

La signification de tout cela peut se traduire en termes profanes ainsi :

1. nous diffusons la structure dans la direction de la structure ($D_1$, premier ordre),
2. nous diffusons la structure dans la direction de la texture ($D_2$, deuxième ordre),
3. nous diffusons la texture dans la direction de la structure ($D_3$, troisième ordre),
4. nous diffusons la texture dans la direction de la texture ($D_4$, quatrième ordre).

Tout coefficient $p$ mis à 0 annule la diffusion, tout coefficient $a$ mis à 0 annule l’anisotropie.

La mise à jour EDP à l’échelle $s$ est donc :

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

où :

• $\kappa$ est le facteur de discrétisation, soit $\frac14$ pour les différences finies centrées,
• $\nu_s$ est le paramètre de régularisation que nous verrons à la section suivante,
• $w_s$ est la pondération d’échelle définie à la section précédente.

Et la resynthèse finale est simplement :

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

Ainsi, si nous résumons l’ensemble de l’algorithme, étant donné $u$ l’image initiale, $u’$ l’image finale, $n$ le nombre final d’échelles :

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

Quelques remarques :

• $n$ n’est pas un paramètre utilisateur mais est déterminé par rapport au $\sigma$ final cible demandé par les paramètres de rayon de l’utilisateur. Cela s’ajuste au niveau de zoom comme sous-produit.
• Le flou gaussien est lui-même la solution 2D isotrope de l’équation de la chaleur : c’est déjà de la diffusion.
• Augmenter le rayon de flou équivaut à laisser la diffusion se dérouler plus longtemps : le signal s’étale plus loin.
• Pour $s > 0$, $H_s$ devient en réalité une différence de gaussiennes. À un certain coefficient de correction d’échelle, la différence de gaussiennes est une approximation du laplacien d’une gaussienne, qui est lui-même une estimation du laplacien à un coefficient d’échelle $\sigma$.
• Appliquer à nouveau un laplacien (éventuellement anisotrope) sur $H_s$ équivaut à une dérivée partielle du 4e ordre (bi-laplacien).

Structure du facteur de correction

Le terme qui contrôle le transport de l’EDP est

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

Dans cette section, nous définirons le facteur $\nu_s$.

Définir des filtres fréquentiels bien comportés pour les photographies

Les photographies sont des reproductions numériques d’une image latente à travers un appareil d’exposition (ouverture du diaphragme, sensibilité ISO du capteur, vitesse d’obturation) et un appareil de discrétisation (ou échantillonnage spatial) (matrice de filtres colorés, grille de pixels). Ce sont des artefacts de la technologie utilisée pour capturer l’image latente et ne concernent pas l’image réelle.

Malheureusement, les heuristiques de capture concernant l’exposition et l’échantillonnage affectent la manière dont nous traitons l’image numérique. J’ai montré sur des disques à la section pochoirs de diffusion 3×3 comment les bords diagonaux se comportent différemment des bords alignés sur la grille (verticaux/horizontaux), même si j’ai choisi les noyaux les plus invariants par rotation : la rotation du contenu de l’image (par rapport à la grille de pixels) va changer la manière dont les gradients discrets sont évalués. L’implication concrète ici est : faire pivoter l’image avant ou après diffuse or sharpen ne produira pas le même résultat.

Mais cela ne s’arrête pas là : l’exposition change aussi la variance du signal. Étant donné un signal blanc $X$, sa variance locale $V_1$ sur une fenêtre d’échantillonnage $\mathcal{N}$ s’exprime :

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

Si, au lieu de capturer $X$, nous avions surexposé la même image d’un facteur $l$, alors la variance de $lX$ deviendrait $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} $$

Ainsi la variance du signal augmente avec le carré du facteur d’exposition. Cela nous importe ici pour deux raisons, qui peuvent être résumées dès maintenant par ceci : tout ce que nous faisons ici consiste à changer la variance de l’image échelle par échelle.

Premièrement, notre étape de floutage par B-spline est une moyenne locale pondérée, et $H_s$ évalué au pixel de coordonnées $(x, y)$ peut en fait s’écrire :

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

avec $h_{i, j}$ les coefficients du noyau B-spline 2D à 5 échantillons, et $X$ le signal flouté de l’échelle précédente (ou l’image initiale pour la première étape). On peut montrer de manière similaire que $H_s$ dépend linéairement du facteur d’exposition. L’équation ci-dessus montre comment $H_s$ peut être vu comme une modulation autour d’une moyenne locale : l’amplitude de cette modulation n’est pas indépendante de l’amplitude du signal. Cela signifie que tout lissage de $H_s$ (réalisé comme un processus diffusif) aura un poids différent et un impact différent sur les détails selon que l’image est surexposée ou sous-exposée, même si le contenu est le même.

Dit autrement, lisser (ou à l’inverse, augmenter la netteté) puis sous-exposer, ou sous-exposer puis lisser n’auront pas le même effet sur les détails, même si l’amplitude globale finale (moyenne) du signal sera la même. Ce n’est pas ce que l’on attend d’un filtre d’image bien comporté : la représentation en données du contenu ne devrait pas affecter la manière dont nous traitons le contenu lui-même. Dans le cadre diffusif, ce n’est pas si dommageable, mais dans le cadre de l’accentuation de la netteté, les détails dans les ombres se retrouvent réellement suraccentués par rapport aux détails dans les hautes lumières, sans aucune forme de normalisation.

Deuxièmement, le floutage par B-spline (ou sa meilleure approximation gaussienne) va changer aussi la variance du signal. Si l’on exprime le signal discret $X$ comme une modulation locale autour de la moyenne globale $\mu$, on obtient $X_n = \mu + \epsilon_n$. Alors, le floutage par B-spline appliqué à $X$ devient :

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

Parce que les coefficients du noyau $h_{k}$ sont normalisés et que $\mu$ est, par définition, constant sur la fenêtre de longueur $k$, on obtient $\mu \sum_{k} h_{k} = \mu$, ce qui signifie que le floutage ne change pas la valeur moyenne. La variance s’exprime alors :

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

À partir de là, on peut montrer que, pour des signaux blancs, non corrélés, $\operatorname{Var}(G_0) = \operatorname{Var}(X) \sum_k h_k^2$, et plus généralement :

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

où $\sum_k h_k^2 = (35 / 128)^2$ pour le filtre B-spline cardinal 2D à 5 échantillons.

Et là encore, nous avons un problème : le même objet échantillonné à une certaine résolution, ou à 4 fois cette résolution, apparaîtrait à la même « fréquence » une étape de décomposition en ondelettes plus tard, ce qui rendrait sa variance de surface $(35 / 128)^2$ fois plus faible. Mais ceci est purement un artefact d’échantillonnage. La variance de l’objet lui-même peut être conceptualisée, en dehors de l’image, dans un cadre continu, comme les modulations locales de couleur autour de sa couleur de surface moyenne. Bien que cette variance idéalisée ne puisse être récupérée par aucun appareil d’imagerie, la manière dont nous traitons le signal devrait au moins être stable en variance de sorte que nous puissions poser que la variance de l’image représente la variance idéalisée de l’objet à un facteur d’échelle constant près.

Cela peut sembler des préoccupations philosophiques jusqu’à ce que nous nous heurtions à un problème pratique de tout logiciel d’imagerie : que se passe-t-il lorsque l’on prévisualise l’effet en zoom avant/arrière ? Comment mettre l’effet à l’échelle de sorte que la prévisualisation réduite reste fidèle au résultat en pleine résolution ?

Ainsi, toutes ces disparités d’échantillonnage doivent être normalisées pour obtenir un filtre d’image qui cherche à manipuler le contenu indépendamment de sa représentation en données.

Définir une métrique de régularisation

Nous avons vu ci-dessus comment la variance du signal est une métrique pertinente pour ce que nous faisons ici : nous pouvons la suivre à travers les étapes de floutage, la relier à l’amplitude du signal, et elle représente la modulation du signal autour de la valeur moyenne.

Malheureusement, nous n’avons pas accès à une métrique de variance une fois que nous entrons dans le schéma de décomposition en ondelettes. Cependant, nous avons vu ci-dessus que $H_s$ était assez proche, conceptuellement, du terme $(\bar{X} - X_i)$ de la variance :

• au lieu d’une moyenne arithmétique, nous utilisons une moyenne pondérée à l’aide des coefficients de B-spline,
• au lieu d’une moyenne globale, nous en utilisons une locale,
• la nature radiale de la B-spline la rend plus invariante par rotation que n’importe quelle moyenne par patch carré.

Nous utiliserons donc l’énergie de la bande $H_s$, évaluée aux mêmes coordonnées de pixel que le pochoir laplacien, définie comme :

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

Pour un signal blanc, non corrélé, à variation lente, $\overline{Q_s} = Q_s / |\mathcal{N}_{3\times 3}|$ devient proche de la variance échelle par échelle et patch par patch.

La régularisation est destinée au problème d’accentuation de la netteté, qui est mal défini : dans ce cadre, nous augmentons l’énergie de chaque couche $H_s$ et il nous faut un paramètre pour la maintenir sous contrôle à un moment donné. C’est une procédure courante dans les problèmes inverses tels que le débruitage et le défloutage, pour lesquels la variation totale est utilisée comme schéma de régularisation depuis un certain temps.

Le modèle de régularisation que nous utiliserons est :

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

avec les paramètres utilisateur :

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

Nous avons montré ci-dessus comment la variance du signal varie avec le carré de la mise à l’échelle de l’exposition, et comment $H_s$ varie linéairement avec la mise à l’échelle de l’exposition. $G_s$ porte la même dépendance linéaire par sa propriété d’être une moyenne locale pondérée.

Ainsi, le rapport $H_s / G_s$ est invariant par exposition. Par identification, $L_s(q) = G_s(q)$ dans l’équation de régularisation, nous utilisons donc l’énergie de bande invariante par exposition :

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

et sa moyenne locale :

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

Normaliser l’échelle et la couverture spatiale

Le pochoir laplacien 3×3 s’étend d’un pas de $2^s$ à chaque échelle $s$, tout comme le fait le noyau B-spline : c’est la base du schéma « à-trous ». L’espace physique couvert par ce noyau augmente avec les échelles d’ondelettes.

Nous pouvons montrer qu’un filtre gaussien 2D de paramètre $\sigma$ change la variance du signal original comme :

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

Étant donné que $4 \pi \sigma^2$ est l’aire effective du disque couvert par le filtre gaussien, ceci correspond à une intuition très simple : le filtre gaussien étale la modulation originale (exprimée sous forme de variance) sur une surface plus grande. C’est la diffusion en un mot. Nous en dérivons :

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

Entre 2 étapes de floutage gaussien, le paramètre de variance du filtre gaussien équivalent varie de :

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

Et d’une échelle à l’autre, ce rayon croît de $\Delta \sigma = \sigma_B 2^s \approx 2^s$, en rappelant que $\sigma_B \approx 1.05\dots$.

Ainsi l’évaluation du laplacien s’étend spatialement au même rythme (même pas à chaque échelle) que les filtres B-spline ou gaussien : le laplacien suit implicitement l’étalement de la variance à travers les échelles, et aucune normalisation d’échelle supplémentaire n’est nécessaire ici.

L’implémentation initiale de diffuse or sharpen avait un boost en $\sigma^2_s$ appliqué au paramètre de régularisation $\lambda$ : en pratique, cela empêche les échelles grossières d’avoir un quelconque impact visible, et les limite à servir de support pour le travail en haute fréquence.

Une autre tentative a été faite pour normaliser $\lambda$ en tenant compte du fait que $Q_s$ a une énergie de bande décroissante à mesure que $s$ augmente, de sorte qu’une métrique d’énergie invariante par échelle a conduit à :

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

Même combiné avec le facteur $\sigma^2_s$ ci-dessus, cela donnait beaucoup trop de poids aux échelles grossières, rendant difficile l’accentuation de la netteté en hautes fréquences alors que les basses dépassaient déjà gravement (donnant un aspect maculé). Ces deux tentatives ont été abandonnées.

Résultat

Après

Avant

L’image ci-dessus paraîtra trop travaillée à la plupart des gens, mais là n’est pas la question : j’ai eu mon lot d’opérateurs d’accentuation défectueux qui fonctionnaient assez bien tant qu’on ne poussait pas la force au-delà de 2 %. De tels filtres d’image n’ont aucun intérêt. On apprend bien plus sur les algorithmes en regardant comment ils échouent qu’en regardant comment ils réussissent dans leur zone de confort.

Cette image est le piège ultime pour tout algorithme d’accentuation de la netteté :

• le premier plan est plus proche et moins brumeux que l’arrière-plan, il réclame donc d’être suraccentué,
• le premier plan proche est aussi beaucoup plus sombre, donc là encore, il est facile de pousser la netteté jusqu’au ciel,
• nous avons une crête de montagne très contrastée qui réclame de produire des halos autour des bords,
• le disque solaire dans un ciel nuageux sera accentué par la plupart des algorithmes avec des bords sombres,
• la quantité de désembuage nécessaire ferait exploser le bruit (ceci a été pris en 2017 avec un Nikon D90 à 200 ISO, nous sommes loin des capteurs actuels).

Ici, nous effectuons donc un débruitage et un défloutage conjoints à grands rayons. Le défloutage utilise une contre-diffusion isotrope (qui décrit le mieux le voile atmosphérique), sur les hautes fréquences comme sur les basses fréquences. Le débruitage utilise une diffusion isophote sur les hautes fréquences, suivant le gradient échantillonné en haute fréquence. Tout cela a été fait sans aucun masquage, dans une seule instance de diffuse or sharpen.

Dans l’ensemble, nous ne voyons aucun dépassement de bord et aucun halo. Le premier plan sombre a été essentiellement ignoré, comme il se doit, mais nous avons récupéré les détails en bas dans la vallée. Nous n’avons obtenu aucun décalage de couleur ni aberration chromatique. Je ne prétendrai pas que cette image est exempte d’artefacts, car, en regardant de près, nous obtenons des traînées de couleur et de nouveaux détails qui pourraient être une véritable reconstruction ou de simples hallucinations du modèle diffusif. Mais il reste que ces artefacts, s’il y en a, ont un aspect suffisamment organique pour passer inaperçus si l’on ne dispose pas de l’image originale.

Cela valide la pertinence du schéma diffusif multi-échelles, ainsi que la stratégie de régularisation.

Après

Avant

Après

Avant

Celle-ci a été prise avec un Sony Alpha 7R Mk 2 mais un Konika Hexanon 40 mm f/1.8. C’est un objectif pancake des années 1970, et la photo est donc 45 Mpx de superbe flou. Les portraits sont moins indulgents que les paysages parce que la peau doit rester saine lorsqu’on récupère la netteté. Ici, nous avons aussi tout le halo d’objectif et les rayons de soleil qui peuvent rapidement dégénérer.

Le détail est particulièrement impressionnant ici parce que nous récupérons les sourcils et les cils sans suraccentuer les cheveux à contre-jour. Dans l’ensemble, il n’y a ni ringing, ni frange, ni halo. Malheureusement, la qualité de la peau s’est dégradée et c’est aussi loin que nous pouvons aller sans masquer la peau pour l’exclure. À partir de là, il s’agit seulement de ramener le contraste global de manière sélective avec une courbe tonale :

Après

Avant

Ce qui est spécifique à Ansel

Les références scientifiques expliquent les briques de base :

• diffusion anisotrope et inpainting par transfert de chaleur,²
• espace d’échelle et diffusion discrète,³⁴
• laplaciens isotropes à 9 points,⁵⁶
• approximation gaussienne par B-spline.⁷

Ce qui est spécifique à Ansel, c’est la manière dont elles sont assemblées :

• une pyramide B-spline à-trous en pleine résolution,
• quatre opérateurs de diffusion paramétrés indépendamment, répartis entre basses et hautes fréquences,
• un régulariseur d’énergie de bande HF en plus normalisé par l’énergie BF locale,
• une sélection d’échelle tenant compte du zoom dans la prévisualisation en chambre noire,
• des mathématiques identiques sur CPU et OpenCL.¹

Le module doit donc être compris comme une synthèse d’ingénierie de plusieurs idées numériques, et non comme l’implémentation littérale d’un seul article.

Perspectives

Le schéma multi-échelles en ondelettes d’EDP de diffusion anisotrope avec régularisation produit des résultats photographiques exploitables, au-delà de la simple preuve de concept. Il permet de tirer parti de l’accentuation et du débruitage conjoints, ainsi que de la diffusion orientée classique. Il permet aussi de rajeunir, par logiciel, de vieux objectifs qui étaient jugés inadaptés à la photographie numérique haute définition. Le module lui-même fournit un terrain de jeu générique d’EDP qui peut être utilisé de nombreuses manières différentes.

Le problème, cependant, est que la nature des réglages est ancrée dans des mathématiques de niveau licence (au moins) et est cryptique pour la plupart des photographes. Expliquer l’effet des paramètres est difficile sans plonger dans ce qu’ils signifient mathématiquement. Essayer de les renommer par leur fonction plutôt que par leur nature est voué à l’échec parce que leur fonction dépend de la manière dont ils sont combinés entre eux, et selon qu’ils sont utilisés dans la plage de valeurs positives ou négatives.

La configuration de diffusion est assez simple et ne nécessite pas de régularisation. Elle n’a pas besoin des 4 ordres tous ensemble, mais seuls les réglages du premier et du troisième ordre suffisent. Dans le cadre isotrope, elle est entièrement équivalente à un floutage gaussien, qui sera moins coûteux à calculer (parce que non itératif).

La configuration d’accentuation, avec débruitage conjoint, est plus compliquée. La seule manière de la rendre plus conviviale est de l’entraîner comme un algorithme d’apprentissage automatique :

• photographier des paires d’images propres et floues/bruitées/brumeuses de la même scène (flou de mouvement, flou de mise au point de l’objectif, objectifs doux),
• effectuer un balayage par force brute des paramètres du module diffuse or sharpen réalisant la reconstruction des images sales, et enregistrer la norme $L_2$ de l’erreur entre les images de référence propres et les reconstructions tentées,
• classer manuellement les images en catégories binaires (floue, bruitée, brumeuse) ou par intensité (le bruit peut être mesuré en PSNR, RMS, … ; le rayon devra probablement être là aussi),
• le problème d’apprentissage automatique devient : dans l’espace 14D des paramètres d’entrée de diffuse or sharpen, quelles sont les 4 directions principales (associées au débruitage, au défloutage, au désembuage, au rayon) qui minimisent l’erreur $E = ||\text{clean} - \text{reconstructed}||_2$ ? Nous cherchons les vecteurs propres de ces directions principales.
• résoudre cela par PLS pondérée (moindres carrés partiels) : $Y = X’ B + c$ pour $n$ paires d’images propres/reconstruites, où :
  • $Y$ est le vecteur $4 × n$ des catégories pour chaque échantillon,
  • $X’$ est le vecteur $14 × n$ standardisé (standardisé composante par composante : $X_j’ = \frac{X_r - \mu}{\sigma}$) des paramètres de diffuse or sharpen,
  • $B$ est la matrice $14 × 4$ mappant nos 14 paramètres cryptiques de D or S vers 4 paramètres conviviaux (c’est l’inconnue ici),
  • $c$ est le résidu (scalaire ou vecteur selon ce qui convient),
  • la norme $L_2$ de l’erreur est utilisée comme pondération de la PLS (probablement injectée dans une exponentielle),
• une fois la matrice $B$ connue (mappant 14D -> 4D), l’inverser pour obtenir le modèle 4D -> 14D,
• ajouter un mode d’interface alternatif exposant les 4 paramètres conviviaux, et une couche interface <-> paramètres convertissant ces 4 vers les 14 arguments d’entrée de D or S (c’est-à-dire : écrire le produit matriciel).

Toute autre tentative de « simplifier » diffuse or sharpen ne sera qu’un stupide travail de réétiquetage obscurcissant le vrai sens des paramètres et empêchant quiconque possédant le bagage mathématique adéquat de le comprendre. Ce serait dommage de faire en sorte que les seules rares personnes capables de le comprendre soient en fait dissuadées ne serait-ce que d’essayer. À l’heure actuelle, D or S est difficile à comprendre, mais au moins il peut être expliqué. Réétiqueter les commandes ne le rendra pas plus facile à comprendre, mais ne fera qu’ajouter une couche supplémentaire de traduction sémantique entre les mathématiques et l’interface, très probablement inexacte et trompeuse de toute façon, ce qui ne fera que le rendre plus exigeant cognitivement à expliquer et à saisir.