Las matemáticas de difuminar o enfocar

Este artículo documenta las matemáticas actualmente implementadas en el módulo diffuse or sharpen de Ansel, tal como se encuentran en src/iop/diffuse.c, src/common/bspline.h, data/kernels/diffuse.cl y data/kernels/bspline.cl. No es una guía de uso. Es una reconstrucción del modelo numérico a partir del código fuente, con las afirmaciones científicas rastreadas hasta las referencias citadas en los comentarios del código.¹

Resumen

El módulo opera sobre una descomposición multiescala à-trous a resolución completa de la imagen, construida a partir de desenfoques de B-spline cardinal. En cada escala, Ansel calcula cuatro operadores anisótropos de segundo orden sobre las bandas de baja frecuencia y de alta frecuencia por separado, regulariza su acción mediante una energía local de la banda de alta frecuencia normalizada por la escala, opcionalmente restringe la actualización a una máscara binaria de reconstrucción (inpainting), y reconstruye la imagen de lo grueso a lo fino.² Los operadores espaciales discretos son plantillas (stencils) de diferencias centradas de 3x3 muestreadas sobre la subred dispersa à-trous de la escala actual, de modo que el mayor soporte físico proviene del paso $2^s$ de la escalera de wavelets y no de un núcleo de EDP mayor.³¹

Motivación

El objetivo del módulo de difusión era simular un pigmento tipo acuarela escapando de una imagen hacia sus bordes, así :

Resultó que simplemente cambiar el signo de la actualización de la ecuación en derivadas parciales de difusión podía en realidad enfocar la imagen. Había trabajado en deconvolución ciega durante años cuando empecé a trabajar en un modelo de difusión, y nunca logré hacer funcionar ni siquiera el método de última generación fuera de los casos bonitos de póster (ISO baja, saturación baja, desenfoque uniforme). Pero esto, por supuesto, necesitaba un parámetro de regularización para evitar la divergencia.

Así fue como nació el módulo diffuse or sharpen, como un marco genérico para simular toda clase de fenómenos difusivos o contra-difusivos.

Modelo continuo

En esencia, el módulo implementa una familia de ecuaciones de difusión anisótropa de la forma

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

donde $u$ es la imagen y $\mathbf{A}$ es un tensor de difusión simétrico definido positivo que dirige el suavizado ya sea a lo largo de las isofotas, a lo largo del gradiente, o de forma isótropa.²⁴

La inspiración original en el código es el modelo de reconstrucción (inpainting) de Qin et al., quienes acoplan la restauración de la estructura y la textura de la imagen mediante transferencia de calor anisótropa.² Ansel conserva ese espíritu, pero aplica la EDP en un dominio de wavelets, por separado sobre las bandas de baja frecuencia y de alta frecuencia, y expone cuatro coeficientes de transporte controlados por el usuario.

Gradientes discretos y anisotropía

Para cada píxel y cada canal, el módulo extrae un vecindario de 3x3 y evalúa las primeras derivadas centradas

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

Este es el gradiente de diferencias centradas estándar usado en los métodos discretos de espacio-escala y de difusión.³

La orientación local del gradiente es entonces :

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

con el habitual valor de reserva de $(1,0)$ cuando la magnitud del gradiente se anula.¹ Nótese que evitaremos calcular el costoso ángulo $\theta$ como la función $\arctan2$ de las componentes $(x, y)$, ya que nunca se usará directamente en lo que sigue.

La intensidad de la anisotropía se convierte del parámetro de usuario $a$ en un coeficiente positivo

$$ \alpha = a^2, $$

mientras que el signo de $a$ selecciona el modo:

• $a = 0$: difusión isótropa,
• $a > 0$: difusión alineada con la dirección de la isofota,
• $a < 0$: difusión alineada con la dirección del gradiente.¹

El término de amortiguamiento es entonces

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

que coincide con la difusión anisótropa de Qin et al. para la reconstrucción (inpainting).²

Para la difusión alineada con la isofota, el tensor escrito en la base local de la imagen es :

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

Para la difusión alineada con el gradiente, Ansel usa la forma invertida :

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

Las plantillas de difusión de 3x3

Dependiendo de los parámetros del usuario, evaluaremos para cada píxel la plantilla $\mathbf{K(A)}$ usando uno de los tensores $\mathbf{A}$ anteriores.

Dado el tensor simétrico $\mathbf{A}$ :

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

construimos el núcleo laplaciano anisótropo, rotado y discreto $\nabla \cdot (\mathbf{A}\nabla u)$ como la plantilla de 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}, $$

salvo por la convención de signos mencionada en el código fuente para los términos fuera de la diagonal.²⁴¹

En el caso isótropo se expande a una constante :

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

que es el laplaciano isótropo clásico al estilo de Oono-Puri. Su interés está en el comportamiento rotacional: comparado con el laplaciano de 5 puntos, el error angular se reduce, lo cual importa cuando la difusión no debe privilegiar los ejes de la imagen.⁵⁶ Este error angular fue evaluado frente a núcleos similares en Rotation-invariant Laplacian for 2D grids .

Para entender el efecto de la dirección en la difusión, empecemos con un disco ruidoso y difuminémoslo con Ansel :

Vale la pena mencionar aquí un par de cosas :

• Solo el desenfoque isótropo produce un disco acromático; los otros retienen algo de ruido de croma a gran escala, porque el ruido de croma crea variaciones locales del gradiente.
• La difusión anisótropa es perfecta cerca de las posiciones norte/este/sur/oeste del círculo, es decir, cuando el ángulo del gradiente está perfectamente alineado con la rejilla de píxeles (0° o 90°). En medio, vemos discrepancias en cuanto a cómo se trata la dirección, debido a las limitaciones numéricas de una plantilla cuadrada de 3×3 rotada.

La pirámide de B-spline à-trous

El análisis multiescala se construye a partir del filtro separable de B-spline cardinal de 5 coeficientes

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

Este filtro es una aproximación gaussiana compacta, razón por la cual aparece de forma natural en los métodos de espacio-escala basados en splines.⁷ Reutilizamos aquí el marco introducido por Johannes Hanika para el módulo contrast equalizer.⁸

La implementación actual usa la escalera à-trous no diezmada histórica:

1. El primer paso bajo a resolución completa es

   $$ G_0 = h_0 * u, $$

2. La banda de detalle más fina es

   $$ H_0 = u - G_0. $$

3. Los niveles más gruesos conservan la misma resolución de imagen y solo agrandan el paso del desenfoque en $2^s$:

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

   lo que significa que en $s=1$, los coeficientes $h_{2}$ se evalúan cada píxel siguiente, en $s=2$, los coeficientes $h_4$ se evalúan cada 4 píxeles, y así sucesivamente…

4. La banda almacenada en la escala $s$ es la diferencia entre dos pasos bajos sucesivos:

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

Así, cada banda vive en la rejilla original de la imagen (sin diezmado), y resolvemos la EDP de lo fino a lo grueso en paralelo, acumulando la solución escala por escala en el búfer de salida hasta que añadimos el residuo final.¹ Este esquema evita apilar la corrección de las escalas finas sobre las escalas gruesas, como es habitual con las pirámides gaussianas o los solucionadores multi-rejilla, y ha mostrado empíricamente una mejor estabilidad en el problema inverso de la reconstrucción de nitidez. Esto equivale a gestionar la energía de cada banda por separado, de un modo que se asemeja a los ecualizadores de audio HiFi, pero dentro de un marco espacial 2D.

Escala del núcleo, escala de banda y la envolvente de escala de la interfaz gráfica

Para estudiar las propiedades del esquema de wavelets, usaremos la propiedad del B-spline cardinal de ser una aproximación de filtro gaussiano. El parámetro gaussiano $\sigma$ controla cómo se componen los núcleos de desenfoque:

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

de modo que las varianzas $\sigma^2$ del núcleo gaussiano se suman bajo convolución.³⁷ En cambio, la varianza de una señal filtrada $X$ :

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

depende del espectro de $X$ y generalmente no es igual a $\sigma^2$. Para desambiguar ambas varianzas, en el resto de este artículo, $\sigma$ se refiere únicamente al parámetro del núcleo gaussiano, nunca a la raíz cuadrada de la varianza de la señal.

Para evitar ambigüedades, también ayuda separar los niveles de paso bajo de las bandas de los niveles:

• $G_0 = h_1 * u$ es el primer paso bajo a resolución completa sobre la imagen de entrada $u$,
• $G_s = h_{2^s} * G_{s-1}$ para $s \ge 1$,
• $H_0 = u - G_0$,
• $H_s = G_{s-1} - G_s$ para $s \ge 1$.

Por lo tanto:

• el índice de paso bajo $s$ se refiere al nivel de desenfoque $G_s$,
• el índice de banda $s$ se refiere a la banda de detalle $H_s$ situada entre $G_{s-1}$ y $G_s$,
• cada nivel permanece muestreado sobre la rejilla original de la imagen.¹

El filtro B-spline aproxima mejor el núcleo gaussiano equivalente de parámetro :

$$ \sigma_B \approx 1.055365. $$

Dado que las varianzas del núcleo gaussiano se suman bajo convolución, el radio de desenfoque efectivo de $G_s$ se deduce directamente de la escalera de análisis à-trous. Por lo tanto, la secuencia acumulativa de paso bajo es

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

con

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

y, en general,

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

Equivalentemente,

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

El código usa $\sigma_{G_s}$ para decidir cuántas escalas se necesitan para igualar el parámetro radius solicitado por el usuario y para ponderar cada banda alrededor del radio central seleccionado por el usuario:

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

donde $z$ es el nivel de zoom del cuarto oscuro (comparado con el raw a resolución completa), $r_c$ es radius_center en la interfaz gráfica, y $r_w$ es radius.¹ Al ajustar estos valores, los usuarios definen el decaimiento de una especie de filtro de “ancho de banda” discretizado centrado en una frecuencia arbitraria. En ese sentido, $r_w$ es el ancho de una envolvente gaussiana en el espacio-escala. No es el radio de la plantilla de la EDP, ni el radio de un único desenfoque. Es la extensión del perfil de ganancia aplicado a las bandas discretas de wavelets.

El esquema de ponderación es estable sin importar el nivel de zoom : los radios se toman en el espacio de la imagen raw (resolución completa). Al previsualizar imágenes reducidas en el cuarto oscuro, las frecuencias más altas quedan truncadas por la reducción de escala, de modo que solo aplicamos la descomposición de wavelets partiendo de la frecuencia más alta disponible y las reponderamos según su radio equivalente a resolución completa. Esto permite una previsualización reducida bastante precisa, aunque pueda ocultar artefactos ruidosos que aparecen en los niveles más altos.

Las bandas cuyo radio de desenfoque equivalente $\sigma_{G_s}$ se sitúa cerca de $r_c / z$ reciben la mayor ganancia, mientras que las bandas más alejadas se atenúan progresivamente. Un $r_w$ pequeño da una selección estrecha de escalas; un $r_w$ grande da una respuesta más amplia y plana sobre las bandas vecinas.

Así que el par (radius_center, radius) debe leerse como un centro y un ancho en el espacio de escalas. El módulo no define directamente una gaussiana en frecuencia de Fourier; en su lugar, define una envolvente gaussiana sobre las bandas multiescala disponibles.¹

Actualización genérica de EDP multiescala

Sean:

• $H_s$ la banda de detalle almacenada en el índice de banda $s$,
• $G_s$ la reconstrucción actual de baja frecuencia usada al resolver la banda $H_s$ durante la síntesis.

En una banda dada $s$, el código construye cuatro respuestas de difusión:

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

donde:

• $K_{1,s}$ y $K_{2,s}$ son plantillas (stencils) de difusión anisótropa 3x3 aplicadas a la reconstrucción actual de baja frecuencia $G_s$,
• $K_{3,s}$ y $K_{4,s}$ son las plantillas análogas aplicadas a la banda de detalle $H_s$,
• $a_1$ a $a_4$ son 4 coeficientes de amortiguación de anisotropía definidos por el usuario (véase más arriba),
• $p_1$ a $p_4$ son 4 coeficientes de actualización (transporte) de la EDP definidos por el usuario (véase más abajo),
• la convolución es por canal sobre la subred dispersa à-trous de paso $2^s$.¹

Todos los núcleos laplacianos $K$ se aplican con la misma distancia de paso que el desenfoque B-spline que se usó para producir la escala de wavelet $s$ sobre la que se aplican, es decir $2^s$. Dado que el núcleo de desenfoque es 5×5 y los núcleos laplacianos son 3×3, eso significa que el laplaciano cubre una cuarta parte de la superficie del desenfoque.

Todos esos $K$ son un esfuerzo por hacer de diffuse or sharpen un marco genérico de EDP multiescala, ya que cruzarán información entre $G_s$ y $H_s$ :

Los gradientes evaluados en $G_s$ tienen más probabilidad de apuntar hacia detalles legítimos de la imagen, y pueden usarse en un contexto de enfoque para ignorar el ruido. Los gradientes evaluados en $H_s$ son más sensibles al ruido y pueden usarse en un contexto difusivo para suavizar el ruido. Separamos la capa sobre la que inferimos la estructura de la imagen (dirección del gradiente) de la capa donde calculamos la actualización de la EDP (laplaciano) y por tanto la difusión.

A falta de un término mejor, esas respuestas $K_i$ están vinculadas a parámetros de la interfaz llamados «órdenes», del primero al cuarto :

• speed, que es el coeficiente de actualización de la EDP,
• anisotropy, que es el coeficiente de amortiguación $a$ del tensor de anisotropía anterior.

Tenemos por tanto 4 parámetros speed, uno para cada orden :

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

Y de forma similar, 4 parámetros anisotropy : $$ (a_1, a_2, a_3, a_4) = (\texttt{first}, \texttt{second}, \texttt{third}, \texttt{fourth}) $$

El significado de todo esto puede traducirse en términos llanos así :

1. difundimos estructura en la dirección de la estructura ($D_1$, primer orden),
2. difundimos estructura en la dirección de la textura ($D_2$, segundo orden),
3. difundimos textura en la dirección de la estructura ($D_3$, tercer orden),
4. difundimos textura en la dirección de la textura ($D_4$, cuarto orden).

Cualquier coeficiente $p$ fijado a 0 cancela la difusión, cualquier coeficiente $a$ fijado a 0 cancela la anisotropía.

La actualización de la EDP en la escala $s$ es por tanto :

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

donde :

• $\kappa$ es el factor de discretización, es decir $\frac14$ para diferencias finitas centradas,
• $\nu_s$ es el parámetro de regularización que veremos en la siguiente sección,
• $w_s$ es la ponderación de escala definida en la sección anterior.

Y la resíntesis final es simplemente :

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

Así pues, si resumimos todo el algoritmo, dados $u$ la imagen inicial, $u’$ la final, $n$ el número final de escalas :

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

Algunas notas:

• $n$ no es un parámetro del usuario, sino que se determina en relación con el $\sigma$ final objetivo solicitado por los parámetros de radio del usuario. Esto se ajusta al nivel de zoom como efecto secundario.
• El desenfoque gaussiano es en sí mismo la solución isótropa 2D de la ecuación del calor : ya es difusión.
• Aumentar el radio de desenfoque equivale a dejar que la difusión se ejecute durante más tiempo : la señal se propaga más lejos.
• Para $s > 0$, $H_s$ se convierte en realidad en una diferencia de gaussianas. Con cierto coeficiente de corrección de escala, la diferencia de gaussianas es una aproximación del laplaciano de una gaussiana, que a su vez es una estimación del laplaciano con un coeficiente de escala $\sigma$.
• Aplicar de nuevo un laplaciano (posiblemente anisótropo) sobre $H_s$ equivale a una derivada parcial de cuarto orden (bilaplaciano).

Estructura del factor de corrección

El término que controla el transporte de la EDP es

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

En esta sección, definiremos el factor $\nu_s$.

Definir filtros de frecuencia bien comportados para fotografías

Las fotografías son reproducciones digitales de una imagen latente a través de un aparato de exposición (apertura del diafragma, sensibilidad ISO del sensor, velocidad de obturación) y un aparato de discretización (o muestreo espacial) (matriz de filtros de color, rejilla de píxeles). Estos son artefactos de la tecnología usada para capturar la imagen latente y no conciernen a la imagen real.

Por desgracia, las heurísticas de captura relativas a la exposición y el muestreo afectan a cómo procesamos la imagen digital. He mostrado en discos, en la sección plantillas de difusión 3×3, cómo los bordes diagonales se comportan de forma distinta a los bordes alineados con la rejilla (verticales/horizontales), aunque elegí los núcleos más invariantes a la rotación : la rotación del contenido de la imagen (comparada con la rejilla de píxeles) cambiará cómo se evalúan los gradientes discretos. La implicación concreta aquí es : rotar la imagen antes o después de diffuse or sharpen no dará el mismo resultado.

Pero no acaba ahí : la exposición cambia también la varianza de la señal. Dada una señal blanca $X$, su varianza local $V_1$ sobre una ventana de muestreo $\mathcal{N}$ se expresa :

$$ \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, en lugar de capturar $X$, sobreexpusiéramos la misma imagen por un factor $l$, entonces la varianza de $lX$ se convertiría en $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} $$

Así que la varianza de la señal aumenta con el cuadrado del factor de exposición. Esto nos importa aquí por dos razones, que pueden resumirse ahora mismo así : todo lo que hacemos aquí es cambiar la varianza por escala de la imagen.

Primero, nuestro paso de desenfoque B-spline es un promedio ponderado local, y $H_s$ evaluado en el píxel de coordenadas $(x, y)$ puede escribirse en realidad :

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

con $h_{i, j}$ los coeficientes del núcleo B-spline 2D de 5 taps, y $X$ la señal desenfocada de la escala anterior (o la imagen inicial en el primer paso). Podemos mostrar de forma similar que $H_s$ depende linealmente del factor de exposición. La ecuación anterior muestra cómo $H_s$ puede verse como una modulación alrededor de un promedio local : la magnitud de esa modulación no es independiente de la magnitud de la señal. Significa que cualquier suavizado de $H_s$ (llevado a cabo como un proceso difusivo) tendrá un peso distinto y un impacto distinto sobre los detalles según si la imagen está sobreexpuesta o subexpuesta, aunque el contenido sea el mismo.

Dicho de otro modo, suavizar (o al contrario, enfocar) y luego subexponer, o subexponer y luego suavizar no tendrá el mismo efecto sobre los detalles, aunque la magnitud global final (el promedio) de la señal sea la misma. Esto no es lo que esperamos de un filtro de imagen bien comportado : la representación en datos del contenido no debería afectar a cómo procesamos el contenido en sí. En el contexto difusivo, no es tan dañino, pero en el contexto de enfoque, los detalles en las sombras quedan realmente sobreenfocados en comparación con los detalles en las luces, sin ningún tipo de normalización.

Segundo, el desenfoque B-spline (o su mejor aproximación gaussiana) cambiará también la varianza de la señal. Si expresamos la señal discreta $X$ como una modulación local alrededor del promedio global $\mu$, obtenemos $X_n = \mu + \epsilon_n$. Entonces, el desenfoque B-spline aplicado a $X$ se convierte en :

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

Como los coeficientes del núcleo $h_{k}$ están normalizados y $\mu$ es, por definición, constante sobre la ventana de longitud $k$, obtenemos $\mu \sum_{k} h_{k} = \mu$, lo que significa que el desenfoque no cambia el valor promedio. La varianza se expresa entonces :

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

A partir de ahí, podemos mostrar que, para señales blancas no correlacionadas, $\operatorname{Var}(G_0) = \operatorname{Var}(X) \sum_k h_k^2$, y de forma más general :

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

donde $\sum_k h_k^2 = (35 / 128)^2$ para el filtro B-spline cardinal 2D de 5 taps.

Y aquí de nuevo tenemos un problema : el mismo objeto muestreado a cierta resolución, o a 4 veces esa resolución, aparecería a la misma «frecuencia» un paso de descomposición en wavelets después, lo que haría que su varianza de superficie fuera $(35 / 128)^2$ menor. Pero esto es puramente un artefacto de muestreo. La varianza del objeto en sí puede conceptualizarse, fuera de la imagen, en un contexto continuo, como las modulaciones locales de color alrededor de su color de superficie promedio. Aunque esta varianza idealizada no puede recuperarse con ningún aparato de imagen, la forma en que manejamos la señal debería al menos ser estable en varianza para poder postular que la varianza de la imagen representa la varianza idealizada del objeto salvo por algún factor de escala constante.

Esto puede parecer una preocupación filosófica hasta que topamos con un problema práctico de cualquier software de imagen : ¿qué ocurre cuando previsualizamos el efecto con zoom de acercamiento/alejamiento ? ¿Cómo escalamos el efecto para que la previsualización reducida siga siendo fiel al resultado a plena resolución ?

Así pues, todas estas discrepancias de muestreo necesitan normalizarse para lograr un filtro de imagen que intente manipular el contenido independientemente de su representación en datos.

Definir una métrica de regularización

Hemos visto más arriba cómo la varianza de la señal es una métrica relevante para lo que hacemos aquí : podemos rastrearla a través de los pasos de desenfoque, vincularla a la magnitud de la señal, y representa la modulación de la señal alrededor del valor promedio.

Por desgracia, no tenemos acceso a una métrica de varianza una vez que entramos en el esquema de descomposición en wavelets. Sin embargo, hemos visto más arriba que $H_s$ era bastante cercano, conceptualmente, al término $(\bar{X} - X_i)$ de la varianza :

• en lugar de una media aritmética, usamos un promedio ponderado con coeficientes B-spline,
• en lugar de un promedio global, usamos uno local,
• la naturaleza radial del B-spline lo hace más invariante a la rotación que cualquier promedio cuadrado por parches.

Así que usaremos la energía de la banda $H_s$, evaluada en las mismas coordenadas de píxel que la plantilla laplaciana, definida como :

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

Para una señal blanca, no correlacionada y de variación lenta, $\overline{Q_s} = Q_s / |\mathcal{N}_{3\times 3}|$ se aproxima a la varianza por escala y por parche.

La regularización está pensada para el problema de enfoque, que está mal definido : en este contexto, aumentamos la energía de cada capa $H_s$ y necesitamos un parámetro para mantenerla a raya en algún punto. Este es un procedimiento habitual en problemas inversos como la reducción de ruido y el desenfoque inverso (deblurring), para los cuales la Variación Total se ha usado como esquema de regularización desde hace bastante tiempo.

El modelo de regularización que usaremos es :

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

con los parámetros del usuario :

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

Hemos mostrado más arriba cómo la varianza de la señal varía con el cuadrado del escalado de exposición, y cómo $H_s$ varía linealmente con el escalado de exposición. $G_s$ arrastra la misma dependencia lineal a través de su propiedad de ser un promedio ponderado local.

Así que la razón $H_s / G_s$ es invariante a la exposición. Por identificación, $L_s(q) = G_s(q)$ en la ecuación de regularización, por lo que usamos la energía de banda invariante a la exposición :

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

y su promedio local :

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

Normalizar la escala y la cobertura espacial

La plantilla laplaciana 3×3 se expande con un paso $2^s$ para cada escala $s$, igual que lo hace el núcleo B-spline : esta es la base del esquema «à-trous». El espacio físico cubierto por ese núcleo aumenta con las escalas de wavelets.

Podemos mostrar que un filtro gaussiano 2D de parámetro $\sigma$ cambia la varianza de la señal original así :

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

Dado que $4 \pi \sigma^2$ es el área efectiva del disco cubierto por el filtro gaussiano, esto coincide con una intuición muy sencilla : el filtro gaussiano propaga la modulación original (expresada como varianza) sobre una superficie mayor. Eso es la difusión en pocas palabras. De ello derivamos :

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

Entre 2 pasos de desenfoque gaussiano, el parámetro de varianza del filtro gaussiano equivalente varía en :

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

Y de una escala a otra, este radio crece en $\Delta \sigma = \sigma_B 2^s \approx 2^s$, recordando que $\sigma_B \approx 1.05\dots$.

Así que la evaluación del laplaciano se expande espacialmente al mismo ritmo (mismo paso en cada escala) que los filtros B-spline o gaussiano : el laplaciano sigue implícitamente la propagación de la varianza a través de las escalas, y aquí no se necesita ninguna normalización de escala adicional.

La implementación inicial de diffuse or sharpen tenía un refuerzo $\sigma^2_s$ aplicado al parámetro de regularización $\lambda$ : en la práctica, eso impide que las escalas gruesas tengan cualquier impacto visible, y las limita a ser el soporte para el trabajo de alta frecuencia.

Se hizo otro intento de normalizar $\lambda$ teniendo en cuenta el hecho de que $Q_s$ tiene una energía de banda decreciente a medida que $s$ aumenta, de modo que una métrica de energía invariante a la escala condujo a :

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

Incluso combinado con el factor $\sigma^2_s$ anterior, eso daba demasiado peso a las escalas gruesas, dificultando el enfoque en altas frecuencias mientras las bajas ya se disparaban gravemente (dando un aspecto emborronado). Ambos intentos se han abandonado.

Resultado

After

Before

La imagen de arriba parecerá recocida a la mayoría, pero eso no viene al caso : he tenido mi buena ración de operadores de enfoque defectuosos que funcionaban bastante bien mientras no subieras la intensidad por encima del 2 %. Tales filtros de imagen no tienen ningún interés. Se aprende mucho más sobre los algoritmos observando cómo fallan que observando cómo triunfan en su punto óptimo.

Esta imagen es la trampa definitiva para cualquier algoritmo de enfoque:

• el primer plano está más cerca y menos brumoso que el fondo, así que pide a gritos que se lo sobreenfoque,
• el primer plano cercano también es mucho más oscuro, así que, de nuevo, es fácil forzar el enfoque hacia el cielo ahí,
• tenemos una cresta montañosa muy contrastada que pide a gritos producir halos alrededor de los bordes,
• el disco solar en un cielo nublado será enfocado por la mayoría de los algoritmos con bordes oscuros,
• la cantidad de eliminación de bruma necesaria haría explotar el ruido (esta foto se tomó en 2017 con una Nikon D90 a 200 ISO, estamos lejos de los sensores actuales).

Así que aquí realizamos una reducción de ruido y desenfoque conjuntos a radios grandes. El desenfoque usa contradifusión isotrópica (que describe mejor la bruma atmosférica), tanto en alta como en baja frecuencia. La reducción de ruido usa difusión de isofotas en alta frecuencia, siguiendo el gradiente muestreado en alta frecuencia. Todo esto se hizo sin enmascaramiento, en una única instancia de diffuse or sharpen.

En general, no vemos sobreoscilación de bordes ni halos. El primer plano oscuro fue en su mayor parte ignorado, como debe ser, pero recuperamos los detalles del valle. No obtuvimos ningún desplazamiento de color ni aberraciones cromáticas. No afirmaré que esta imagen esté libre de artefactos, porque, al mirar de cerca, obtenemos vetas de color y nuevos detalles que podrían ser una reconstrucción real o simples alucinaciones del modelo difusivo. Pero el hecho sigue siendo que estos artefactos, si los hay, se ven lo bastante orgánicos como para pasar desapercibidos si no tenemos la imagen original.

Esto valida la pertinencia del esquema difusivo multiescala, junto con la estrategia de regularización.

After

Before

After

Before

Esta se tomó con una Sony Alpha 7R Mk 2 pero con un Konika Hexanon 40 mm f/1.8. Es un objetivo pancake de los años 1970, y por tanto la imagen son 45 Mpx de hermoso desenfoque. Los retratos son menos indulgentes que los paisajes porque la piel necesita mantenerse sana al recuperar el enfoque. Aquí también tenemos toda la reflexión del objetivo y los rayos de sol que pueden degenerar rápidamente.

El detalle es especialmente impresionante aquí porque recuperamos las cejas y las pestañas sin sobreenfocar el cabello a contraluz. En general, no hay ringing, franjas de color ni halos. Por desgracia, la calidad de la piel se ha degradado y es lo más lejos que podemos llegar sin enmascarar la piel para excluirla. A partir de ahí, solo se trata de recuperar el contraste global de forma selectiva con una curva de tonos:

After

Before

Qué es específico de Ansel

Las referencias científicas explican los bloques de construcción:

• difusión anisotrópica e inpainting por transferencia de calor,²
• espacio-escala y difusión discreta,³⁴
• laplacianos isotrópicos de 9 puntos,⁵⁶
• aproximación gaussiana por B-splines.⁷

Lo específico de Ansel es la forma en que se ensamblan:

• una pirámide B-spline à-trous a resolución completa,
• cuatro operadores de difusión parametrizados independientemente, repartidos entre bajas y altas frecuencias,
• un regularizador de energía de banda de HF adicionalmente normalizado por la energía LF local,
• selección de escala consciente del zoom en la vista previa del cuarto oscuro,
• matemáticas idénticas en CPU y OpenCL.¹

Así que el módulo debe entenderse como una síntesis de ingeniería de varias ideas numéricas, no como una implementación literal de un único artículo.

Perspectivas

El esquema de wavelets multiescala de PDE de difusión anisotrópica con regularización produce resultados fotográficos utilizables, más allá de la mera prueba de concepto. Permite aprovechar el enfoque y la reducción de ruido conjuntos, junto con la difusión orientada habitual. También permite rejuvenecer, mediante software, objetivos antiguos que se consideraban inadecuados para la fotografía digital de alta definición. El módulo en sí proporciona un campo de pruebas de PDE genérico que puede usarse de muchas maneras diferentes.

El problema, sin embargo, es que la naturaleza de los ajustes se fundamenta en matemáticas de (al menos) nivel de grado y resulta críptica para la mayoría de los fotógrafos. Explicar el efecto de los parámetros es difícil sin sumergirse en lo que significan matemáticamente. Intentar renombrarlos por su función en lugar de por su naturaleza está condenado al fracaso, porque su función depende de cómo se combinan entre sí, y de si se usan en el rango de valores positivos o negativos.

La configuración de difusión es bastante sencilla y no requiere regularización. No necesita los 4 órdenes en conjunto, sino que basta con los ajustes de primer y tercer orden. En la configuración isotrópica, es totalmente equivalente a un desenfoque gaussiano, que será menos costoso de calcular (por ser no iterativo).

La configuración de enfoque, junto con la reducción de ruido conjunta, es más complicada. La única forma de hacerla más fácil de usar es entrenarla como un algoritmo de aprendizaje automático:

• fotografiar pares de imágenes limpias y borrosas/ruidosas/brumosas de la misma escena (desenfoque de movimiento, desenfoque de desenfoque del objetivo, objetivos blandos),
• realizar un barrido de parámetros por fuerza bruta del módulo diffuse or sharpen que reconstruya las imágenes sucias, y registrar la norma $L_2$ del error entre las imágenes de referencia limpias y las reconstrucciones intentadas,
• clasificar manualmente las imágenes en categorías binarias (borrosa, ruidosa, brumosa) o por intensidad (el ruido puede medirse como PSNR, RMS, …; probablemente el radio también deba estar ahí),
• el problema de aprendizaje automático se convierte en: en el espacio 14D de parámetros de entrada de diffuse or sharpen, ¿cuáles son las 4 direcciones principales (asociadas a reducción de ruido, desenfoque, eliminación de bruma, radio) que minimizan el error $E = ||\text{clean} - \text{reconstructed}||_2$? Buscamos los vectores propios de esas direcciones principales.
• resolverlo mediante PLS ponderado (mínimos cuadrados parciales): $Y = X’ B + c$ para $n$ pares de imágenes limpias/reconstruidas, donde:
  • $Y$ es el vector $4 × n$ de categorías para cada muestra,
  • $X’$ es el vector $14 × n$ estandarizado (estandarizado componente a componente: $X_j’ = \frac{X_r - \mu}{\sigma}$) de parámetros de diffuse or sharpen,
  • $B$ es la matriz $14 × 4$ que mapea nuestros 14 crípticos parámetros de D or S a 4 parámetros fáciles de usar (es la incógnita aquí),
  • $c$ es el residuo (escalar o vector según lo que encaje),
  • la norma $L_2$ del error se usa como ponderación del PLS (probablemente inyectada en una exponencial),
• una vez conocida la matriz $B$ (que mapea 14D -> 4D), invertirla para obtener el modelo 4D -> 14D,
• añadir un modo de GUI alternativo que exponga los 4 parámetros fáciles de usar, y una capa GUI <-> parámetros que convierta esos 4 en los 14 argumentos de entrada de D or S (es decir: escribir el producto matricial).

Cualquier otro intento de “simplificar” diffuse or sharpen no será más que un trabajo tonto de reetiquetado que ofusca el significado real de los parámetros e impide que cualquiera con la formación matemática adecuada lo entienda. Sería una lástima asegurarse de que las pocas personas capaces de entenderlo acaben disuadidas de siquiera intentarlo. Ahora mismo, D or S es difícil de entender, pero al menos puede explicarse. Reetiquetar los controles no lo hará más fácil de entender, sino que solo añadirá una capa extra de traducción semántica entre las matemáticas y la GUI, muy probablemente inexacta y engañosa de todos modos, lo que solo hará que sea más exigente a nivel cognitivo explicarlo y comprenderlo.