A matemática do difundir ou aguçar

Este artigo documenta a matemática atualmente implementada no módulo diffuse or sharpen do Ansel, tal como se encontra em src/iop/diffuse.c, src/common/bspline.h, data/kernels/diffuse.cl e data/kernels/bspline.cl. Não é um guia de uso. É uma reconstrução do modelo numérico a partir do código-fonte, com as afirmações científicas rastreadas até as referências citadas nos comentários do código.¹

Resumo

O módulo opera sobre uma decomposição multiescala à-trous em resolução plena da imagem, construída a partir de borramentos por B-spline cardinal. Em cada escala, o Ansel calcula quatro operadores anisotrópicos de segunda ordem separadamente sobre as bandas de baixa frequência e de alta frequência, regulariza sua ação por uma energia local da banda de alta frequência normalizada pela escala, opcionalmente restringe a atualização a uma máscara binária de inpainting e reconstrói a imagem do grosseiro ao fino.² Os operadores espaciais discretos são estênceis 3x3 de diferenças centrais amostrados na sub-rede esparsa à-trous da escala atual, de modo que o suporte físico maior vem do passo $2^s$ da escada de wavelets, e não de um núcleo de EDP maior.³¹

Motivação

O objetivo do módulo de difusão era simular um pigmento semelhante a aquarela escapando de uma imagem para suas bordas, assim :

Acabou se revelando que simplesmente mudar o sinal da atualização da equação diferencial parcial de difusão podia, na verdade, aguçar a imagem. Eu havia trabalhado por anos em deconvolução cega quando comecei a trabalhar em um modelo de difusão, e nunca consegui fazer funcionar nem mesmo o método estado da arte fora dos belos casos de pôster (ISO baixo, saturação baixa, borramento uniforme). Mas isso, é claro, precisava de um parâmetro de regularização para evitar a divergência.

Foi assim que nasceu o módulo difundir ou aguçar, como um arcabouço genérico para a simulação de todo tipo de fenômeno difusivo ou contradifusivo.

Modelo contínuo

No cerne, o módulo implementa uma família de equações de difusão anisotrópica da forma

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

onde $u$ é a imagem e $\mathbf{A}$ é um tensor de difusão simétrico positivo que direciona a suavização seja ao longo das isófotas, ao longo do gradiente ou isotropicamente.²⁴

A inspiração original no código é o modelo de inpainting de Qin et al., que acopla a restauração da estrutura e da textura da imagem por meio de transferência de calor anisotrópica.² O Ansel mantém esse espírito, mas aplica a EDP em um domínio de wavelets, separadamente sobre as bandas de baixa frequência e de alta frequência, e expõe quatro coeficientes de transporte controlados pelo usuário.

Gradientes discretos e anisotropia

Para cada pixel e cada canal, o módulo extrai uma vizinhança 3x3 e avalia primeiras 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 é o gradiente por diferenças centrais padrão usado nos métodos discretos de espaço-escala e de difusão.³

A orientação local do gradiente é então :

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

com o recurso usual de retornar a $(1,0)$ quando a magnitude do gradiente se anula.¹ Note que evitaremos calcular o dispendioso ângulo $\theta$ como a função $\arctan2$ das componentes $(x, y)$, já que ele nunca será usado diretamente no que se segue.

A intensidade da anisotropia é convertida do parâmetro do usuário $a$ em um coeficiente positivo

$$ \alpha = a^2, $$

enquanto o sinal de $a$ seleciona o modo:

• $a = 0$: difusão isotrópica,
• $a > 0$: difusão alinhada com a direção da isófota,
• $a < 0$: difusão alinhada com a direção do gradiente.¹

O termo de amortecimento é então

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

que corresponde à difusão anisotrópica de Qin et al. para inpainting.²

Para a difusão alinhada com a isófota, o tensor escrito na base local da imagem é :

$$ \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 a difusão alinhada com o gradiente, o Ansel usa a 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}. $$

Os estênceis de difusão 3x3

Dependendo dos parâmetros do usuário, avaliaremos para cada pixel o estêncil $\mathbf{K(A)}$ usando um dos tensores $\mathbf{A}$ acima.

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

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

construímos o núcleo laplaciano anisotrópico, rotacionado e discreto $\nabla \cdot (\mathbf{A}\nabla u)$ como o estêncil 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}, $$

a menos da convenção de sinal mencionada no código-fonte para os termos fora da diagonal.²⁴¹

No caso isotrópico, expande-se para uma constante :

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

que é o laplaciano isotrópico clássico ao estilo Oono-Puri. Seu interesse está no comportamento rotacional: comparado ao laplaciano de 5 pontos, o erro angular é reduzido, o que importa quando a difusão não deve privilegiar os eixos da imagem.⁵⁶ Esse erro angular foi avaliado em comparação a núcleos similares em Rotation-invariant Laplacian for 2D grids .

Para entender o efeito da direção na difusão, vamos começar com um disco ruidoso e difundi-lo com o Ansel :

Vale a pena mencionar algumas coisas aqui :

• Apenas o borramento isótropo produz um disco acromático; os outros retêm algum ruído de croma em larga escala, porque o ruído de croma cria variações locais de gradiente.
• A difusão anisótropa é perfeita perto das posições norte/leste/sul/oeste do círculo, isto é, quando o ângulo do gradiente está perfeitamente alinhado à grade de pixels (0° ou 90°). No meio-termo, vemos discrepâncias quanto à forma como a direção é tratada, devidas às limitações numéricas de um estêncil quadrado 3×3 rotacionado.

A pirâmide de B-spline à-trous

A análise multiescala é construída a partir do filtro B-spline cardinal separável de 5 tomadas

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

Este filtro é uma aproximação gaussiana compacta, e é por isso que ele aparece naturalmente nos métodos de espaço-escala baseados em splines.⁷ Reutilizamos aqui o arcabouço introduzido por Johannes Hanika para o módulo contrast equalizer.⁸

A implementação atual usa a escada à-trous não decimada histórica:

1. O primeiro passa-baixa em resolução plena é

   $$ G_0 = h_0 * u, $$

2. A banda de detalhe mais fina é

   $$ H_0 = u - G_0. $$

3. Os níveis mais grosseiros mantêm a mesma resolução da imagem e apenas ampliam o passo do borramento por $2^s$:

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

   significando que em $s=1$, as tomadas $h_{2}$ são avaliadas a cada pixel seguinte, em $s=2$, as tomadas $h_4$ são avaliadas a cada 4 pixels, e assim por diante…

4. A banda armazenada na escala $s$ é a diferença entre dois passa-baixa sucessivos:

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

Assim, cada banda vive na grade original da imagem (sem decimação), e resolvemos a EDP do fino ao grosseiro em paralelo, acumulando a solução escala por escala no buffer de saída até adicionarmos o resíduo final.¹ Este esquema evita empilhar a correção das escalas finas sobre as escalas grosseiras, como é costume com pirâmides gaussianas ou solucionadores multigrade, e mostrou empiricamente melhor estabilidade no problema inverso da reconstrução de nitidez. Isso é equivalente a gerenciar a energia de cada banda separadamente, de um modo que se assemelha aos equalizadores de áudio HiFi, mas dentro de um arcabouço espacial 2D.

Escala do núcleo, escala da banda e o envelope de escala da GUI

Para estudar as propriedades do esquema de wavelets, usaremos a propriedade da B-spline cardinal de ser uma aproximação de filtro gaussiano. O parâmetro gaussiano $\sigma$ controla como os núcleos de borramento se compõem:

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

de modo que as variâncias $\sigma^2$ do núcleo gaussiano se somam sob convolução.³⁷ Em contraste, a variância de um sinal filtrado $X$ :

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

depende do espectro de $X$ e geralmente não é igual a $\sigma^2$. Para desambiguar ambas as variâncias, no restante deste artigo, $\sigma$ refere-se apenas ao parâmetro do núcleo gaussiano, nunca à raiz quadrada da variância do sinal.

Para evitar ambiguidade, também ajuda separar os níveis de passa-baixa das bandas dos níveis:

• $G_0 = h_1 * u$ é o primeiro passa-baixa em resolução plena sobre a imagem 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$.

Assim:

• o índice de passa-baixa $s$ refere-se ao nível de borramento $G_s$,
• o índice de banda $s$ refere-se à banda de detalhe $H_s$ situada entre $G_{s-1}$ e $G_s$,
• cada nível permanece amostrado na grade original da imagem.¹

O filtro B-spline aproxima melhor o núcleo gaussiano equivalente de parâmetro :

$$ \sigma_B \approx 1.055365. $$

Como as variâncias do núcleo gaussiano se somam sob convolução, o raio de borramento efetivo de $G_s$ segue diretamente da escada de análise à-trous. Portanto, a sequência cumulativa de passa-baixa é

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

com

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

e, em geral,

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

O código usa $\sigma_{G_s}$ para decidir quantas escalas são necessárias para corresponder ao parâmetro radius solicitado pelo usuário e para ponderar cada banda em torno do raio central selecionado pelo usuário:

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

onde $z$ é o nível de zoom da sala escura (comparado ao raw em resolução plena), $r_c$ é radius_center na GUI e $r_w$ é radius.¹ Ao definir esses valores, os usuários determinam o decaimento de uma espécie de filtro de “largura de banda” discretizado, centrado em uma frequência arbitrária. Nesse sentido, $r_w$ é a largura de um envelope gaussiano no espaço-escala. Não é o raio do estêncil da EDP, nem o raio de um único borramento. É a dispersão do perfil de ganho aplicado às bandas discretas de wavelets.

O esquema de ponderação é estável independentemente do nível de zoom : os raios são tomados no espaço da imagem raw (resolução plena). Ao pré-visualizar imagens reduzidas na sala escura, as frequências mais altas são truncadas pela redução de escala, então aplicamos a decomposição de wavelets apenas a partir da frequência mais alta disponível e as reponderamos de acordo com seu raio equivalente em resolução plena. Isso permite uma pré-visualização reduzida bastante precisa, ainda que possa ocultar artefatos ruidosos que aparecem nos níveis mais altos.

As bandas cujo raio de desfoque equivalente $\sigma_{G_s}$ está próximo de $r_c / z$ recebem o maior ganho, enquanto as bandas mais distantes são progressivamente atenuadas. Um $r_w$ pequeno dá uma seleção estreita de escalas; um $r_w$ grande dá uma resposta mais ampla e mais plana sobre as bandas vizinhas.

Portanto, o par (radius_center, radius) deve ser lido como um centro e uma largura no espaço de escalas. O módulo não define diretamente uma gaussiana na frequência de Fourier; em vez disso, ele define uma envoltória gaussiana sobre as bandas multiescala disponíveis.¹

Atualização genérica multiescala da EDP

Seja:

• $H_s$ a banda de detalhe armazenada no índice de banda $s$,
• $G_s$ a reconstrução de baixa frequência atual usada ao resolver a banda $H_s$ durante a síntese.

Em uma dada banda $s$, o código constrói quatro respostas de difusão:

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

onde:

• $K_{1,s}$ e $K_{2,s}$ são estênceis de difusão anisotrópica 3x3 aplicados à reconstrução de baixa frequência atual $G_s$,
• $K_{3,s}$ e $K_{4,s}$ são os estênceis análogos aplicados à banda de detalhe $H_s$,
• $a_1$ a $a_4$ são 4 coeficientes de amortecimento de anisotropia definidos pelo usuário (ver acima),
• $p_1$ a $p_4$ são 4 coeficientes de atualização (transporte) da EDP definidos pelo usuário (ver abaixo),
• a convolução é feita por canal sobre a sub-rede esparsa à-trous de passo $2^s$.¹

Todos os núcleos laplacianos $K$ são aplicados com a mesma distância de passo do desfoque B-spline que foi usado para produzir a escala wavelet $s$ sobre a qual são aplicados, ou seja, $2^s$. Dado que o núcleo de desfoque é 5×5 e os núcleos laplacianos são 3×3, isso significa que o laplaciano cobre um quarto da superfície do desfoque.

Todos esses $K$ são um esforço para tornar o diffuse or sharpen um framework genérico de EDP multiescala, pois eles cruzam informação entre $G_s$ e $H_s$ :

Os gradientes avaliados sobre $G_s$ têm mais probabilidade de apontar para detalhes legítimos da imagem, e podem ser usados em um cenário de aumento de nitidez para ignorar o ruído. Os gradientes avaliados sobre $H_s$ são mais sensíveis ao ruído e podem ser usados em um cenário difusivo para suavizar o ruído. Separamos a camada sobre a qual inferimos a estrutura da imagem (direção do gradiente) da camada onde calculamos a atualização da EDP (laplaciano) e, portanto, a difusão.

Por falta de um termo melhor, essas respostas $K_i$ estão ligadas a parâmetros da interface chamados “ordens”, da primeira à quarta :

• speed, que é o coeficiente de atualização da EDP,
• anisotropy, que é o coeficiente de amortecimento $a$ do tensor de anisotropia acima.

Temos, portanto, 4 parâmetros speed, um para cada ordem :

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

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

O significado de tudo isso pode ser traduzido em termos leigos assim :

1. difundimos estrutura na direção da estrutura ($D_1$, primeira ordem),
2. difundimos estrutura na direção da textura ($D_2$, segunda ordem),
3. difundimos textura na direção da estrutura ($D_3$, terceira ordem),
4. difundimos textura na direção da textura ($D_4$, quarta ordem).

Qualquer coeficiente $p$ definido como 0 cancela a difusão, qualquer coeficiente $a$ definido como 0 cancela a anisotropia.

A atualização da EDP na escala $s$ é, portanto :

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

onde :

• $\kappa$ é o fator de discretização, ou seja, $\frac14$ para diferenças finitas centradas,
• $\nu_s$ é o parâmetro de regularização que veremos na próxima seção,
• $w_s$ é a ponderação de escala definida na seção anterior.

E a ressíntese final é simplesmente :

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

Então, se resumirmos todo o algoritmo, dados $u$ a imagem inicial, $u’$ a final, $n$ o 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*} $$

Algumas observações:

• $n$ não é um parâmetro do usuário, mas é determinado em relação ao $\sigma$ final alvo solicitado pelos parâmetros de raio do usuário. Isso se ajusta ao nível de zoom como subproduto.
• O desfoque gaussiano é, ele próprio, a solução isótropa 2D da equação do calor : ele já é difusão.
• Aumentar o raio de desfoque equivale a deixar a difusão correr por mais tempo : o sinal se espalha mais longe.
• Para $s > 0$, $H_s$ na verdade se torna uma diferença de gaussianas. A partir de algum coeficiente de correção de escala, a diferença de gaussianas é uma aproximação do laplaciano de uma gaussiana, que por sua vez é uma estimativa do laplaciano a um coeficiente de escala $\sigma$.
• Aplicar novamente um laplaciano (possivelmente anisótropo) sobre $H_s$ equivale a uma derivada parcial de 4ª ordem (bi-laplaciano).

Estrutura do fator de correção

O termo que controla o transporte da EDP é

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

Nesta seção, definiremos o fator $\nu_s$.

Definindo filtros de frequência bem comportados para fotografias

Fotografias são reproduções digitais de uma imagem latente através de um aparato de exposição (abertura do diafragma, sensibilidade ISO do sensor, velocidade do obturador) e de um aparato de discretização (ou amostragem espacial) (matriz de filtros de cor, grade de pixels). Esses são artefatos da tecnologia usada para capturar a imagem latente e não dizem respeito à imagem em si.

Infelizmente, as heurísticas de captura relativas à exposição e à amostragem afetam a forma como processamos a imagem digital. Mostrei, com discos, na seção estênceis de difusão 3×3, como as arestas diagonais se comportam de forma diferente das arestas alinhadas à grade (verticais/horizontais), mesmo tendo eu escolhido os núcleos mais invariantes a rotação : a rotação do conteúdo da imagem (em relação à grade de pixels) mudará como os gradientes discretos são avaliados. A implicação concreta aqui é : rotacionar a imagem antes ou depois do diffuse or sharpen não produzirá o mesmo resultado.

Mas não para por aí : a exposição também altera a variância do sinal. Dado um sinal branco $X$, sua variância local $V_1$ sobre uma janela de amostragem $\mathcal{N}$ é expressa :

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

Se, em vez de capturar $X$, superexpuséssemos a mesma imagem por um fator $l$, então a variância de $lX$ se tornaria $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} $$

Assim, a variância do sinal aumenta com o quadrado do fator de exposição. Isso nos importa aqui por duas razões, que podem ser resumidas agora mesmo assim : tudo o que fazemos aqui é alterar a variância por escala da imagem.

Primeiro, nossa etapa de desfoque B-spline é uma média local ponderada, e $H_s$ avaliado no pixel de coordenadas $(x, y)$ pode na verdade ser escrito :

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

com $h_{i, j}$ os coeficientes do núcleo B-spline 2D de 5 taps, e $X$ o sinal desfocado da escala anterior (ou a imagem inicial para a primeira etapa). Podemos mostrar de forma semelhante que $H_s$ depende linearmente do fator de exposição. A equação acima mostra como $H_s$ pode ser visto como uma modulação em torno de uma média local : a magnitude dessa modulação não é independente da magnitude do sinal. Isso significa que qualquer suavização de $H_s$ (realizada como um processo difusivo) terá um peso diferente e um impacto diferente sobre os detalhes conforme a imagem esteja superexposta ou subexposta, mesmo que o conteúdo seja o mesmo.

Dito de outra forma, suavizar (ou, ao contrário, aumentar a nitidez) e depois subexpor, ou subexpor e depois suavizar não terá o mesmo efeito sobre os detalhes, mesmo que a magnitude final geral (média) do sinal seja a mesma. Não é isso que esperamos de um filtro de imagem bem comportado : a representação de dados do conteúdo não deveria afetar como processamos o conteúdo em si. No cenário difusivo, isso não é tão prejudicial, mas no cenário de aumento de nitidez, os detalhes nas sombras ficam realmente com nitidez excessiva em comparação com os detalhes nas realces, sem nenhum tipo de normalização.

Segundo, o desfoque B-spline (ou sua melhor aproximação gaussiana) também alterará a variância do sinal. Se expressarmos o sinal discreto $X$ como uma modulação local em torno da média global $\mu$, obtemos $X_n = \mu + \epsilon_n$. Então, o desfoque B-spline aplicado a $X$ torna-se :

$$ \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 os coeficientes do núcleo $h_{k}$ são normalizados e $\mu$ é, por definição, constante sobre a janela de comprimento $k$, obtemos $\mu \sum_{k} h_{k} = \mu$, o que significa que o desfoque não altera o valor médio. A variância é então expressa :

$$ \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 daí, podemos mostrar que, para sinais brancos e não correlacionados, $\operatorname{Var}(G_0) = \operatorname{Var}(X) \sum_k h_k^2$, e mais geralmente :

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

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

E aqui novamente temos um problema : o mesmo objeto amostrado a alguma resolução, ou a 4 vezes essa resolução, apareceria na mesma “frequência” uma etapa de decomposição em wavelets depois, o que tornaria sua variância de superfície $(35 / 128)^2$ menor. Mas isso é puramente um artefato de amostragem. A variância do próprio objeto pode ser conceitualizada, fora da imagem, em um cenário contínuo, como as modulações locais de cor em torno de sua cor média de superfície. Embora essa variância idealizada não possa ser recuperada por nenhum aparato de imageamento, a forma como tratamos o sinal deveria ao menos ser estável em variância, para que possamos postular que a variância da imagem representa a variância idealizada do objeto a menos de algum fator de escala constante.

Isso pode parecer uma preocupação filosófica até esbarrarmos em um problema prático de qualquer software de imageamento : o que acontece quando visualizamos o efeito com zoom ampliado/reduzido ? Como escalonamos o efeito para que a pré-visualização reduzida ainda seja fiel ao resultado em resolução plena ?

Portanto, todas essas discrepâncias de amostragem precisam ser normalizadas para se obter um filtro de imagem que tente manipular o conteúdo independentemente de sua representação de dados.

Definindo uma métrica de regularização

Vimos acima como a variância do sinal é uma métrica relevante para o que estamos fazendo aqui : podemos rastreá-la através das etapas de desfoque, ligá-la à magnitude do sinal, e ela representa a modulação do sinal em torno do valor médio.

Infelizmente, não temos acesso a uma métrica de variância uma vez que entramos no esquema de decomposição em wavelets. No entanto, vimos acima que $H_s$ estava, conceitualmente, bem próximo do termo $(\bar{X} - X_i)$ da variância :

• em vez de uma média aritmética, usamos uma média ponderada com coeficientes B-spline,
• em vez de uma média global, usamos uma local,
• a natureza radial da B-spline a torna mais invariante a rotação do que qualquer média por bloco quadrado.

Então usaremos a energia da banda $H_s$, avaliada nas mesmas coordenadas de pixel que o estêncil laplaciano, definida como :

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

Para um sinal branco, não correlacionado e de variação lenta, $\overline{Q_s} = Q_s / |\mathcal{N}_{3\times 3}|$ torna-se próximo da variância por escala e por bloco.

A regularização destina-se ao problema de aumento de nitidez, que é mal definido : nesse cenário, aumentamos a energia de cada camada $H_s$ e precisamos de um parâmetro para contê-la em algum ponto. Este é um procedimento comum em problemas inversos como redução de ruído e deblurring, para os quais a Variação Total tem sido usada como esquema de regularização há bastante tempo.

O modelo de regularização que usaremos é :

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

com os parâmetros do usuário :

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

Mostramos acima como a variância do sinal varia com o quadrado do escalonamento de exposição, e como $H_s$ varia linearmente com o escalonamento de exposição. $G_s$ carrega a mesma dependência linear por meio de sua propriedade de ser uma média local ponderada.

Portanto, a razão $H_s / G_s$ é invariante à exposição. Por identificação, $L_s(q) = G_s(q)$ na equação de regularização, então usamos a energia de banda invariante à exposição :

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

e sua média 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 $$

Normalizando a escala e a cobertura espacial

O estêncil laplaciano 3×3 se expande por um passo de $2^s$ a cada escala $s$, assim como o núcleo B-spline : essa é a base do esquema “à-trous”. O espaço físico coberto por esse núcleo aumenta com as escalas de wavelets.

Podemos mostrar que um filtro gaussiano 2D de parâmetro $\sigma$ altera a variância do sinal original assim :

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

Dado que $4 \pi \sigma^2$ é a área efetiva do disco coberta pelo filtro gaussiano, isso corresponde a uma intuição muito simples : o filtro gaussiano espalha a modulação original (expressa como variância) sobre uma superfície maior. Isso é difusão em poucas palavras. Derivamos disso :

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

Entre 2 etapas de desfoque gaussiano, o parâmetro de variância do filtro gaussiano equivalente varia de :

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

E de uma escala para outra, esse raio cresce por $\Delta \sigma = \sigma_B 2^s \approx 2^s$, lembrando que $\sigma_B \approx 1.05\dots$.

Portanto, a avaliação do laplaciano se expande espacialmente à mesma taxa (mesmo passo em cada escala) que os filtros B-spline ou gaussianos : o laplaciano segue implicitamente o espalhamento da variância através das escalas, e nenhuma normalização de escala adicional é necessária aqui.

A implementação inicial do diffuse or sharpen tinha um reforço de $\sigma^2_s$ aplicado ao parâmetro de regularização $\lambda$ : na prática, isso impede que as escalas grosseiras tenham qualquer impacto visível, e as limita a servir de suporte para o trabalho em alta frequência.

Outra tentativa foi feita para normalizar $\lambda$ levando em conta o fato de que $Q_s$ tem uma energia de banda decrescente à medida que $s$ aumenta, de modo que uma métrica de energia invariante à escala levou a :

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

Mesmo combinada com o fator $\sigma^2_s$ acima, isso resultava em peso excessivo demais nas escalas grosseiras, dificultando o aumento de nitidez em altas frequências enquanto as baixas já estavam ultrapassando os limites de forma grave (dando uma aparência borrada). Ambas as tentativas foram abandonadas.

Resultado

After

Before

A imagem acima vai parecer exagerada para a maioria, mas isso não vem ao caso: já tive minha cota de operadores de nitidez defeituosos que funcionavam bem o suficiente desde que você não aumentasse a intensidade acima de 2%. Filtros de imagem assim têm interesse zero. Você aprende muito mais sobre algoritmos observando como eles falham do que como eles têm sucesso no seu ponto ideal.

Esta imagem é a armadilha definitiva para todos os algoritmos de nitidez:

• o primeiro plano está mais próximo e menos enevoado que o fundo, então está implorando para ser supernitidificado,
• o primeiro plano próximo também é muito mais escuro, então, de novo, é fácil aumentar demais a nitidez ali no céu,
• temos uma crista de montanha muito contrastada implorando para produzir halos ao redor das bordas,
• o disco do sol num céu nublado será nitidificado pela maioria dos algoritmos com bordas escuras,
• a quantidade de remoção de névoa necessária faria o ruído explodir (esta foi tirada em 2017 com uma Nikon D90 a 200 ISO, estamos longe dos sensores atuais).

Então, aqui realizamos uma redução de ruído e um desborramento conjuntos em raios grandes. O desborramento usa contradifusão isotrópica (que descreve melhor o enevoamento atmosférico), tanto em altas frequências quanto em baixas frequências. A redução de ruído usa difusão de isófota em altas frequências, seguindo o gradiente amostrado em altas frequências. Tudo isso foi feito sem nenhum mascaramento, numa única instância do diffuse or sharpen.

No geral, não vemos nenhum estouro de borda nem halos. O primeiro plano escuro foi majoritariamente ignorado, como deveria, mas recuperamos os detalhes lá embaixo no vale. Não obtivemos desvio de cor nem aberrações cromáticas. Não vou afirmar que esta imagem esteja livre de artefatos, porque, ao olhar de perto, obtemos raios de cor e novos detalhes que podem ser reconstrução real ou puras alucinações do modelo difusivo. Mas o ponto permanece: esses artefatos, se houver, parecem orgânicos o suficiente para passar despercebidos se não tivermos a imagem original.

Isso valida a relevância do esquema difusivo multiescala, juntamente com a estratégia de regularização.

After

Before

After

Before

Esta foi tirada com uma Sony Alpha 7R Mk 2, mas com uma Konika Hexanon 40 mm f/1.8. É uma lente pancake dos anos 1970, e por isso a imagem é 45 Mpx de um belo borrão. Retratos são menos tolerantes que paisagens, porque a pele precisa se manter saudável ao trazer a nitidez de volta. Aqui também temos todo o lensflare e os raios de sol que podem degenerar rapidamente.

O detalhe é especialmente impressionante aqui porque recuperamos as sobrancelhas e os cílios sem supernitidificar o cabelo em contraluz. No geral, não há ringing, fringing nem haloing. Infelizmente, a qualidade da pele degradou e é o máximo que podemos ir sem mascarar a pele para excluí-la. A partir daí, é só trazer de volta o contraste global seletivamente com uma curva de tons:

After

Before

O que é específico do Ansel

As referências científicas explicam os blocos de construção:

• difusão anisotrópica e inpainting por transferência de calor,²
• espaço-escala e difusão discreta,³⁴
• Laplacianos isotrópicos de 9 pontos,⁵⁶
• aproximação Gaussiana por B-spline.⁷

O que é específico do Ansel é a maneira como eles são montados:

• uma pirâmide B-spline à trous em resolução plena,
• quatro operadores de difusão parametrizados independentemente, divididos entre baixas e altas frequências,
• um regularizador de energia de banda de HF adicionalmente normalizado pela energia local de LF,
• seleção de escala ciente do zoom na pré-visualização da sala escura,
• matemática idêntica em CPU e OpenCL.¹

Portanto, o módulo deve ser entendido como uma síntese de engenharia de várias ideias numéricas, e não como uma implementação literal de um único artigo.

Perspectivas

O esquema de wavelets multiescala de EDP de difusão anisotrópica com regularização produz resultados fotográficos utilizáveis, para além da mera prova de conceito. Ele permite aproveitar nitidez e redução de ruído conjuntas, junto com a difusão orientada regular. Também permite rejuvenescer, por software, lentes antigas que eram consideradas inadequadas para a fotografia digital de alta definição. O próprio módulo fornece um playground genérico de EDP que pode ser usado de muitas maneiras diferentes.

O problema, porém, é que a natureza das configurações está fundamentada em matemática (pelo menos) de nível de graduação e é enigmática para a maioria dos fotógrafos. Explicar o efeito dos parâmetros é difícil sem mergulhar no que eles significam matematicamente. Tentar renomeá-los pela sua função em vez de pela sua natureza está fadado ao fracasso, porque a função deles depende de como são combinados entre si e de se são usados na faixa de valores positivos ou negativos.

A configuração de difusão é bastante direta e não requer regularização. Ela não precisa das 4 ordens todas juntas, mas apenas as configurações de primeira e terceira ordem são suficientes. Na configuração isotrópica, é totalmente equivalente a um borrão Gaussiano, que será menos custoso de computar (por ser não iterativo).

A configuração de nitidez, junto com a redução de ruído conjunta, é mais complicada. A única maneira de torná-la mais amigável ao usuário é treiná-la como um algoritmo de aprendizado de máquina:

• fotografar pares de imagens limpas e borradas/ruidosas/enevoadas da mesma cena (borrão de movimento, borrão de desfoque da lente, lentes suaves),
• realizar uma varredura de parâmetros por força bruta do módulo diffuse or sharpen efetuando a reconstrução das imagens sujas, e registrar a norma $L_2$ do erro entre as imagens de referência limpas e as reconstruções tentadas,
• classificar manualmente as imagens em categorias binárias (borrada, ruidosa, enevoada) ou por intensidade (o ruído pode ser medido como PSNR, RMS, …; o raio provavelmente também precisará estar lá),
• o problema de aprendizado de máquina torna-se: no espaço 14D dos parâmetros de entrada do diffuse or sharpen, quais são as 4 direções principais (associadas à redução de ruído, ao desborramento, à remoção de névoa, ao raio) que minimizam o erro $E = ||\text{clean} - \text{reconstructed}||_2$? Estamos procurando os autovetores dessas direções principais.
• resolver isso por PLS ponderado (mínimos quadrados parciais): $Y = X’ B + c$ para $n$ pares de imagens limpas/reconstruídas, onde:
  • $Y$ é o vetor $4 × n$ de categorias para cada amostra,
  • $X’$ é o vetor $14 × n$ padronizado (padronizado componente a componente: $X_j’ = \frac{X_r - \mu}{\sigma}$) dos parâmetros do diffuse or sharpen,
  • $B$ é a matriz $14 × 4$ que mapeia nossos 14 parâmetros enigmáticos do D or S para 4 parâmetros amigáveis ao usuário (é a incógnita aqui),
  • $c$ é o resíduo (escalar ou vetor, dependendo do que melhor se ajustar),
  • a norma $L_2$ do erro é usada como a ponderação do PLS (provavelmente injetada numa exponencial),
• uma vez conhecida a matriz $B$ (mapeando 14D -> 4D), inverta-a para obter o modelo 4D -> 14D,
• adicionar um modo de GUI alternativo expondo os 4 parâmetros amigáveis ao usuário, e uma camada GUI <-> parâmetros convertendo esses 4 nos 14 argumentos de entrada do D or S (ou seja: escrever o produto matricial).

Qualquer outra tentativa de “simplificar” o diffuse or sharpen será apenas um trabalho tolo de reetiquetagem, ofuscando o significado real dos parâmetros e impedindo que qualquer pessoa com a formação matemática adequada o entenda. Seria uma pena garantir que as poucas pessoas capazes de entendê-lo fossem, na verdade, desencorajadas até de tentar. No momento, o D or S é difícil de entender, mas ao menos pode ser explicado. Reetiquetar os controles não o tornará mais fácil de entender, mas apenas adicionará uma camada extra de tradução semântica entre a matemática e a GUI, muito provavelmente imprecisa e enganosa de qualquer forma, o que só tornará mais cognitivamente exigente explicá-lo e compreendê-lo.