Teoria · Forma · Fourier elíptico
Uma semente é uma soma de elipses.
Dê uma volta no contorno de uma semente anotando x e y. As duas listas se repetem a cada volta, e o que se repete vira série de Fourier. Escrita com x e y separados, cada termo da série desenha uma elipse, e meia dúzia de elipses já descreve a forma.
Aprofundar · matemática e referências
O contorno como curva
O contorno de uma semente é uma curva fechada. Para fazer conta com ela, percorra o contorno em velocidade constante a partir de um ponto qualquer e chame de \(t\) o comprimento já andado. Depois de uma volta, \(t\) chega ao perímetro \(T\) e a caminhada recomeça no mesmo lugar.
Definição
Um contorno parametrizado pelo comprimento de arco é um par de funções \(x(t)\) e \(y(t)\), com \(0 \le t \lt T\), tal que o ponto \(\big(x(t), y(t)\big)\) anda sobre o contorno com velocidade 1 e volta ao início em \(t = T\). As duas funções são periódicas: \(x(t+T) = x(t)\) e \(y(t+T) = y(t)\).
Na prática o contorno chega como polígono, a lista de vértices que sai da segmentação. Entre dois vértices, \(x\) e \(y\) variam em linha reta com \(t\). São funções contínuas, com quinas nos vértices, e mais adiante fica claro que isso basta para a série de Fourier convergir para elas em todos os pontos.
Duas séries de Fourier
Cada coordenada tem a sua série. Com \(\theta = 2\pi t/T\), o ângulo que mede a fração da volta já percorrida,
com os coeficientes de sempre:
e \(c_n\), \(d_n\), \(C_0\) iguais, trocando \(x\) por \(y\). O par \((A_0, C_0)\) é o centro do contorno, a média dos pontos ao longo do perímetro. O termo de ordem \(n\) dá \(n\) voltas enquanto o ponto dá uma: \(n = 1\) descreve a forma grossa, os termos altos descrevem os detalhes da borda.
Cada termo é uma elipse
O truque de Kuhl e Giardina [1] é olhar o termo \(n\) das duas séries ao mesmo tempo, como uma matriz 2 × 2 aplicada a um ponto do círculo:
Resultado 1 · o harmônico é uma elipse
Quando \(t\) dá uma volta no contorno, o harmônico \(n\) percorre \(n\) vezes uma elipse centrada na origem. Os semieixos são os valores singulares \(\sigma_1 \ge \sigma_2\) de \(M_n\), a área é \(\pi\,\lvert\det M_n\rvert\) e o sentido do giro é o sinal de \(\det M_n\).
Esboço. O ponto \((\cos n\theta, \sin n\theta)\) anda sobre o círculo unitário. Escreva a decomposição em valores singulares \(M_n = U\Sigma V^\top\): \(V^\top\) gira o círculo, \(\Sigma\) estica os eixos por \(\sigma_1\) e \(\sigma_2\) e \(U\) gira de novo. O círculo vira uma elipse de semieixos \(\sigma_1\) e \(\sigma_2\), e a área é multiplicada por \(\lvert\det M_n\rvert = \sigma_1\sigma_2\). Um determinante negativo inverte o sentido.
A curva inteira é a soma das elipses. O centro da segunda elipse anda sobre a primeira, o da terceira anda sobre a segunda, e a ponta da última desenha o contorno. O instrumento mais abaixo desenha essa cadeia girando sobre uma semente real.
Resultado 2 · a área sai dos coeficientes
A área dentro do contorno, com sinal positivo quando o contorno é percorrido no sentido anti-horário, é
Esboço. Pelo teorema de Green, \(A = \tfrac12\oint (x\,dy - y\,dx)\). Substitua as séries (1) e integre em \(\theta\) de 0 a \(2\pi\). Os produtos de termos de ordens diferentes integram a zero, e cada ordem \(n\) deixa \(\pi n (a_n d_n - b_n c_n)\).
Na semente da Figura 1, a fórmula (4) com 40 harmônicos dá 727,12 px², e a fórmula do cadarço aplicada aos vértices do polígono dá 727,14 px². A diferença, 0,003%, é a cauda da série que ficou de fora.
Do polígono aos coeficientes
As integrais (2) parecem pedir integração numérica. Para um polígono não pedem. Sejam \(K\) lados, com \(\Delta x_p\) e \(\Delta y_p\) os deslocamentos do lado \(p\), \(\Delta t_p = \sqrt{\Delta x_p^2 + \Delta y_p^2}\) o seu comprimento, \(t_p = \Delta t_1 + \dots + \Delta t_p\) e \(\theta_p = 2\pi t_p / T\). Kuhl e Giardina [1] chegaram a
e \(c_n\), \(d_n\) iguais com \(\Delta y_p\) no lugar de \(\Delta x_p\). O centro é a média de cada lado pesada pelo seu comprimento,
que é a forma fechada de (2) para uma função linear em cada lado, equivalente à que está no artigo original.
De onde vem (5)
Integre (2) por partes. Como \(x\) é periódica, o termo de fronteira some e fica \(a_n = -\frac{1}{\pi n}\int_0^T x'(t)\sin n\theta\,dt\). Num lado do polígono, \(x'(t) = \Delta x_p/\Delta t_p\) é constante, e a integral de \(\sin n\theta\) de \(t_{p-1}\) a \(t_p\) vale \(-\frac{T}{2\pi n}\big(\cos n\theta_p - \cos n\theta_{p-1}\big)\). Somando os lados sai (5).
Duas consequências práticas. A conta é exata para o polígono, sem reamostrar e sem integral numérica, e custa \(K \cdot N\) senos e cossenos. E o fator \(1/n^2\) mostra que os coeficientes de um polígono caem pelo menos como \(1/n^2\). A soma dos seus módulos converge, e então a série converge uniformemente [9]: a reconstrução encosta no polígono em todos os pontos, quinas incluídas.
Reconstruir com N harmônicos
Cortar as séries em \(N\) termos dá a curva \(z_N(t) = \big(x_N(t), y_N(t)\big)\), a soma do centro com as \(N\) primeiras elipses. Ela guarda \(4N + 2\) números. A pergunta é quanto se perde.
Resultado 3 · truncar é a melhor aproximação
Entre todas as curvas cujas coordenadas são polinômios trigonométricos de grau \(N\), a série truncada é a mais próxima do contorno no sentido do erro quadrático médio, e esse erro é a cauda dos coeficientes:
Esboço. As funções \(1, \cos n\theta, \sin n\theta\) são ortogonais em \([0, T)\). Cortar a série é projetar \(x\) e \(y\) ortogonalmente no espaço dos polinômios trigonométricos de grau \(N\), e a projeção ortogonal é o ponto mais perto. O que sobra é ortogonal ao que fica, e o teorema de Pitágoras nesse espaço (a identidade de Parseval) dá (7) [9].
A cauda tem infinitos termos, mas não é preciso somá-la. A energia total \(\tfrac{1}{T}\int_0^T \lvert z - z_0 \rvert^2 dt\), com \(z_0 = (A_0, C_0)\), sai exata do polígono, lado por lado, e \(e_N^2\) é essa energia menos metade da soma dos quadrados dos \(N\) primeiros harmônicos. É assim que o instrumento calcula o erro. O mesmo enunciado volta na página da SVD com outro nome: lá, truncar uma decomposição ortogonal dá a melhor aproximação de posto baixo.
Com \(N = 1\) sai a elipse que melhor resume a semente, e o erro é 6,3% do comprimento. Com \(N = 3\) aparece o recorte em forma de rim, e o erro cai para 1,9%. Com 10 harmônicos a curva fica a 0,48% do contorno, e com 20, a 0,21%.
Na mediana das 145 sementes, o erro cai de 3,45% com um harmônico para 1,40% com três e 0,83% com cinco. Metade das sementes fica abaixo de 1% com \(N = 5\), e 90% delas com \(N = 10\). A mais exigente, uma Calligonum mongolicum de contorno recortado, só chega lá com 17 harmônicos, e a Bromus inermis da linha tracejada, com 14.
Resultado 4 · forma simétrica pelo centro não tem harmônico par
Se o contorno é simétrico pelo centro, isto é, \(z(t + T/2) - z_0 = -\big(z(t) - z_0\big)\), então \(a_n = b_n = c_n = d_n = 0\) para todo \(n\) par.
Esboço. Troque \(t\) por \(t + T/2\) em (2). O cosseno ganha o fator \(\cos(n\theta + n\pi) = (-1)^n \cos n\theta\) e a função troca de sinal em torno do centro, logo \(a_n = -(-1)^n a_n\). Para \(n\) par isso obriga \(a_n = 0\). O mesmo vale para \(b_n\), \(c_n\) e \(d_n\).
Grão de arroz é quase simétrico pelo centro. Nas 40 silhuetas de arroz, o segundo harmônico tem, na mediana, 1,35% do comprimento, e o terceiro, 2,97%. Por isso o erro do arroz quase não se mexe de \(N = 1\) para \(N = 2\) (de 3,57% para 3,23%) e despenca com \(N = 3\) (1,40%). É o mesmo degrau que a mediana da Figura 2 mostra entre 2 e 3.
Instrumento · elipses sobre um contorno real
Tamanho de cada harmônico, \(\sqrt{(a_n^2+b_n^2+c_n^2+d_n^2)/2}\), em % do comprimento (escala log)
Os coeficientes são calculados aqui, no seu navegador, pela fórmula (5), a partir dos 128 pontos de cada contorno embutidos na página. Tracejado é o contorno; laranja, a soma das N elipses; em azul, as elipses encadeadas, cada uma centrada na ponta da anterior. Com "normalizar", a figura é girada, escalada e transladada para a primeira elipse ficar deitada, com semieixo maior 1, e o ponto que anda parte da ponta do eixo maior. Na soja, a razão dos semieixos fica perto de 1 e a primeira elipse quase vira círculo.
Normalizar pela primeira elipse
A mesma semente, girada ou com o contorno começando em outro ponto, dá coeficientes diferentes. Para comparar formas é preciso tirar da conta a posição, o tamanho, a rotação e o ponto de partida. Kuhl e Giardina fazem isso usando só a primeira elipse [1]. Primeiro, o ângulo \(\theta_1\) que leva o início do contorno para a ponta do eixo maior da primeira elipse:
Depois, com \(R(\varphi)\) a rotação do plano pelo ângulo \(\varphi\), cada harmônico é deslocado no tempo por \(n\theta_1\), girado pelo ângulo \(\psi_1\) do eixo maior e dividido pelo semieixo maior \(E\):
Resultado 5 · normalizar é uma SVD 2 × 2
Depois de (9), a primeira matriz vira \(M_1'' = \begin{pmatrix} 1 & 0 \\ 0 & \pm\sigma_2/\sigma_1 \end{pmatrix}\): a primeira elipse fica deitada no eixo \(x\), com semieixo maior 1, e o contorno começa na ponta desse eixo.
Esboço. Escreva \(M_1 = U\Sigma V^\top\). O ponto da primeira elipse mais longe do centro corresponde à direção \(v_1\), primeira coluna de \(V\), e (8) é o ângulo de \(v_1\), o autovetor principal de \(M_1^\top M_1\). Então \(M_1 R(\theta_1)\) tem primeira coluna \(\sigma_1 u_1\), \(\psi_1\) é o ângulo de \(u_1\) e \(E = \sigma_1\). Girar por \(-\psi_1\) e dividir por \(\sigma_1\) deixa \(\operatorname{diag}(1, \pm\sigma_2/\sigma_1)\). É a mesma decomposição da página seguinte, em tamanho pequeno.
O que se perde fica contado. Saem a posição \((A_0, C_0)\), o tamanho \(E\), a rotação \(\psi_1\) e o ponto de partida \(\theta_1\), e três números da primeira elipse ficam fixos (\(a_1'' = 1\), \(b_1'' = c_1'' = 0\)). Dos \(4N + 2\) números sobram \(4N - 3\), que descrevem só a forma. Ficam duas ambiguidades. Trocar \(\theta_1\) por \(\theta_1 + \pi\) também leva o início para uma ponta do eixo maior, e as duas escolhas diferem pelo sinal dos harmônicos pares. E a semente que cai com a outra face para cima aparece espelhada, com outros coeficientes: a normalização não reconhece o espelho como a mesma forma.
Conferimos a invariância com a semente da Figura 1. Girando o contorno 40°, ampliando 1,7 vez e começando a contagem a um terço da volta, os coeficientes normalizados mudaram no máximo 6 × 10⁻¹⁶, o tamanho do arredondamento da máquina. Os três primeiros harmônicos normalizados são estes:
| n | \(a_n''\) | \(b_n''\) | \(c_n''\) | \(d_n''\) |
|---|---|---|---|---|
| 1 | 1 | 0 | 0 | 0,490 |
| 2 | 0,024 | −0,049 | −0,138 | 0,036 |
| 3 | 0,086 | 0,013 | 0,089 | 0,026 |
O \(d_1'' = 0{,}490\) é a razão entre os semieixos da primeira elipse. O segundo harmônico, que não existiria numa forma simétrica pelo centro, tem \(c_2'' = -0{,}138\): é a assimetria entre o lado côncavo e o lado convexo do rim.
Onde falha · semente quase redonda
A normalização depende da direção do eixo maior da primeira elipse, e essa direção só está bem definida quando \(\sigma_1\) é bem maior que \(\sigma_2\). Uma perturbação pequena na matriz pode girar os vetores singulares de um ângulo da ordem da perturbação dividida por \(\sigma_1 - \sigma_2\) [10]. Em 24 das 33 silhuetas de soja, o semieixo menor da primeira elipse passa de 90% do maior (mediana 0,928). Trocar o desfoque usado para tirar o contorno de 0,6 para 1,0 px gira a primeira elipse da soja em 1,6° na mediana e em até 9,0°. No arroz, a mesma troca gira 0,05°, e no máximo 0,2°. Na soja, os coeficientes normalizados carregam esse giro. Para sementes quase redondas, compare quantidades que não dependem do ângulo, como os semieixos \(\sigma_1\) e \(\sigma_2\) de cada harmônico, ou normalize por uma referência marcada na semente.
Os descritores clássicos
Antes de 1982 havia duas famílias de descritores de Fourier. Zahn e Roskies [2] expandiram em série o ângulo da tangente ao longo do contorno, uma função que não muda com translação nem com escala e que a rotação só desloca por uma constante. O preço é depender da direção da tangente, sensível à escada de pixels, e a série truncada não garante uma curva fechada. Granlund [3] e Persoon e Fu [4] trataram o contorno como número complexo \(z = x + iy\), com uma série só:
Cada termo é um círculo de raio \(\lvert Z_k \rvert\) que gira \(k\) vezes por volta, no sentido anti-horário para \(k \gt 0\) e no horário para \(k \lt 0\). É o que o vídeo mostra, num contorno-modelo.
Resultado 6 · uma elipse são dois círculos
O harmônico elíptico \(n\) e o par de termos complexos \(\pm n\) guardam a mesma informação:
e os semieixos da elipse são \(\sigma_1 = \lvert Z_n\rvert + \lvert Z_{-n}\rvert\) e \(\sigma_2 = \big\lvert \lvert Z_n\rvert - \lvert Z_{-n}\rvert \big\rvert\).
Esboço. Em (3), escreva \(\cos n\theta = (e^{in\theta} + e^{-in\theta})/2\) e \(\sin n\theta = (e^{in\theta} - e^{-in\theta})/2i\) e agrupe os termos em \(e^{\pm in\theta}\). Para os semieixos, confira que \(\lvert Z_n\rvert^2 - \lvert Z_{-n}\rvert^2 = a_n d_n - b_n c_n = \det M_n\) e que \(2(\lvert Z_n\rvert^2 + \lvert Z_{-n}\rvert^2) = a_n^2 + b_n^2 + c_n^2 + d_n^2 = \sigma_1^2 + \sigma_2^2\).
Os três termos do começo do vídeo (\(k = -1, 0, 1\)) são o centro e uma elipse. Os 25 do final são o centro e 12 elipses. A escrita elíptica separa o que cada eixo faz e dá uma normalização com sentido geométrico, e foi a que se espalhou na morfometria de plantas. Rohlf e Archie [5] compararam os métodos de Fourier em asas de mosquito e acharam os descritores elípticos os mais promissores. O programa SHAPE [6] e o pacote Momocs [7] fazem a conta desta página a partir de imagens. Em sementes de trigo, os descritores elípticos acharam mais regiões do genoma ligadas à forma do grão do que a razão entre comprimento e largura [8].
No SeedCounter
No SeedCounter
A máquina propõe o contorno de cada semente e a pessoa confere. O contorno conferido é um polígono, e o polígono é a entrada exata da fórmula (5): os coeficientes saem dele sem reamostrar e sem integral numérica. O app mede cada semente sobre um contorno simplificado de até 48 lados. Com 48 vértices há 48 valores de \(x\) e 48 de \(y\), e cada harmônico gasta dois de cada. A partir da ordem 24, a série passa a descrever as quinas do polígono, não a semente. O Fourier elíptico ainda não existe no app, nem como cálculo nem como coluna da planilha. Esta página mostra a conta que eles fariam e quantos harmônicos fazem sentido para cada resolução.
Onde falha · Fourier copia o contorno que recebe
A série descreve com fidelidade a forma que entrou, inclusive a errada. Se a sombra colou na borda ou duas sementes encostadas viraram uma, os coeficientes descrevem a sombra e o par. O post A luz decide o contorno mostra de onde vem boa parte desses erros, e por isso a conferência vem antes da conta.
A resolução também manda. Nas silhuetas desta página, com cerca de 51 px de comprimento, a mediana do erro fica abaixo de meio pixel a partir de \(N = 5\). O que os harmônicos seguintes acrescentam é menor que a resolução da imagem e passa a descrever a escada de pixels e o desfoque. E a série compara pontos pelo comprimento de arco, não por pontos homólogos: o hilo de duas sementes pode cair em valores de \(t\) diferentes.
Dados
As 145 silhuetas de 64 × 64 px dos exemplos públicos do SeedCounter (arroz, soja, café, milho e espécies nativas), as mesmas da página Aprender. Cada silhueta vem deitada pelo eixo maior e escalada para cerca de 51 px de comprimento, então as medidas estão em px da silhueta, não em mm. O contorno é tirado da maior região de cada silhueta, com buracos preenchidos, depois de um desfoque gaussiano de 0,6 px, e reamostrado em 128 pontos igualmente espaçados. O comprimento é o diâmetro de Feret máximo. O script que refaz as figuras está em _src/figuras/teoria-forma-fourier.py.
Referências
- Kuhl F.P., Giardina C.R. (1982). Elliptic Fourier features of a closed contour. Computer Graphics and Image Processing 18(3), 236–258. doi:10.1016/0146-664X(82)90034-X
- Zahn C.T., Roskies R.Z. (1972). Fourier descriptors for plane closed curves. IEEE Transactions on Computers C-21(3), 269–281. doi:10.1109/TC.1972.5008949
- Granlund G.H. (1972). Fourier preprocessing for hand print character recognition. IEEE Transactions on Computers C-21(2), 195–201. doi:10.1109/TC.1972.5008926
- Persoon E., Fu K.S. (1977). Shape discrimination using Fourier descriptors. IEEE Transactions on Systems, Man, and Cybernetics 7(3), 170–179. doi:10.1109/TSMC.1977.4309681
- Rohlf F.J., Archie J.W. (1984). A comparison of Fourier methods for the description of wing shape in mosquitoes (Diptera: Culicidae). Systematic Zoology 33(3), 302–317. doi:10.2307/2413076
- Iwata H., Ukai Y. (2002). SHAPE: a computer program package for quantitative evaluation of biological shapes based on elliptic Fourier descriptors. Journal of Heredity 93(5), 384–385. doi:10.1093/jhered/93.5.384
- Bonhomme V., Picq S., Gaucherel C., Claude J. (2014). Momocs: outline analysis using R. Journal of Statistical Software 56(13). doi:10.18637/jss.v056.i13
- Williams K., Munkvold J., Sorrells M. (2013). Comparison of digital image analysis using elliptic Fourier descriptors and major dimensions to phenotype seed shape in hexaploid wheat (Triticum aestivum L.). Euphytica 190(1), 99–116. doi:10.1007/s10681-012-0783-0
- Stein E.M., Shakarchi R. (2003). Fourier Analysis: An Introduction. Princeton University Press.
- Golub G.H., Van Loan C.F. (2013). Matrix Computations, 4ª ed. Johns Hopkins University Press.
Primeiro o contorno, depois as elipses.
O SeedCounter propõe o contorno de cada semente e você confere. A forma só vira número depois disso.