16

Para um trabalho de uma disciplina, eu construi o seguinte código em Python com uso dos pacotes Matplotlib e scikit-image:

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np
from scipy import interpolate
from mpl_toolkits.mplot3d import Axes3D
from pylab import *
from skimage.filter import gabor_kernel
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

# Função utilizada para reescalar o kernel
def resize_kernel(aKernelIn, iNewSize):
    x = np.array([v for v in range(len(aKernelIn))])
    y = np.array([v for v in range(len(aKernelIn))])
    z = aKernelIn

    xx = np.linspace(x.min(), x.max(), iNewSize)
    yy = np.linspace(y.min(), y.max(), iNewSize)

    aKernelOut = np.zeros((iNewSize, iNewSize), np.float)    
    oNewKernel = interpolate.RectBivariateSpline(x, y, z)
    aKernelOut = oNewKernel(xx, yy)

    return aKernelOut


if __name__ == "__main__":
    fLambda = 3.0000000001   # comprimento de onda (em pixels)
    fTheta  = 0              # orientação (em radianos)
    fSigma  = 0.56 * fLambda # envelope gaussiano (com 1 oitava de largura de banda)
    fPsi    = np.pi / 2      # deslocamento (offset)

    # Tamanho do kernel (3 desvios para cada lado, para limitar cut-off)
    iLen = int(math.ceil(3.0 * fSigma))
    if(iLen % 2 != 0):
        iLen += 1

    # Obtém o kernel Gabor para os parâmetros definidos
    z = np.real(gabor_kernel(fLambda, theta=fTheta, sigma_x=fSigma, sigma_y=fSigma, offset=fPsi))

    # Plotagem do kernel

    fig = plt.figure(figsize=(16, 9))
    fig.suptitle(r'Gabor kernel para $\lambda=3.0$, $\theta=0.0$, $\sigma=0.56\lambda$ e $\psi=\frac{\pi}{2}$', fontsize=25)
    grid = gridspec.GridSpec(1, 2, width_ratios=[1, 2])

    # Gráfico 2D
    plt.gray()
    ax = fig.add_subplot(grid[0])
    ax.set_title(u'Visualização 2D')
    ax.imshow(z, interpolation='nearest')
    ax.set_xticklabels([v for v in range((-iLen/2)-1, iLen/2+1)], fontsize=10)
    ax.set_yticklabels([v for v in range((-iLen/2)-1, iLen/2+1)], fontsize=10)

    # Gráfico em 3D

    # Reescalona o kernel para uma exibição melhor
    z = resize_kernel(z, 300)
    # Eixos x e y no intervalo do tamanho do kernel
    x = np.linspace(-iLen/2, iLen/2, 300)
    y = x
    x, y = meshgrid(x, y)

    ax = fig.add_subplot(grid[1], projection='3d')
    ax.set_title(u'Visualização 3D')
    ax.plot_surface(x, y, z, cmap='hot')

    print(iLen)

    plt.show()

Note a definição da variável fLambda na linha 29:

fLambda = 3.0000000001

Se a variável é definida dessa forma (isto é, com um grande número de decimais, mas com valor bem próximo de 0), a plotagem do kernel Gabor (que é a intenção desse programa) aparece como esperado:

inserir a descrição da imagem aqui

Porém, se a variável é definida em uma das formas abaixo (isto é, com o valor redondo '3'):

fLambda = 3.0
# ou
fLambda = float(3)

[...] a plotagem aparece com um resultado bastante diferente do esperado:

inserir a descrição da imagem aqui

A minha primeira impressão era de que a interpolação utilizada no reescalonamento do kernel (função resize_kernel) estava errada, mas pode-se observar na visualização 2D (que utiliza diretamente o Kernel devolvido pelo pacote scikit-image, isto é, SEM interpolação), que o kernel já está diferente.

Outra informação: ao usar o valor 3.000000000000001 (15 decimais) o resultado é igual ao da primeira imagem; ao usar valores com mais decimais (por exemplo, 3.0000000000000001 - 1 casa decimal a mais), o resultado já passa a ser igual à segunda imagem.

Esse resultado parece correto pra alguém (será o Gabor kernel tão diferente para valores inteiros e de ponto flutuante assim)? Pode haver algum erro de precisão envolvido? Será alguma especificidade do Python ou mesmo do pacote scikit-image?

Versão do Python: 2.7.5 (windows 7 32bit) Versão do Matplotlib: 1.3.1 Versão do numpy: 1.8.0 Versão do scikit-image: 0.9.3

15
  • 3
    Eu "apostaria" num erro de precisão. Lembre-se que 3.0000000001 é maior que 3, se em qualquer momento for exigido um teto ele será 4 no primeiro caso e 3 no segundo (quando vi o ceil no seu código quase achei que havia achado o erro, mas aí percebi que ele envolvia fSigma e não fLambda). E, claro, ao aumentar demais o limite de casas decimais você atinge o limite da representação em ponto flutuante, fazendo o valor voltar a ser 3 (e, portanto, o resultado ser igual ao da segunda imagem).
    – mgibsonbr
    Commented 14/12/2013 às 11:52
  • Obriado pelo comentário, @mgibsonbr. Você deve ter razão, talvez o pacote scikit-image use um ceil internamente e cause diferenças no resultado. Estou fazendo um teste semelhante em C++ com OpenCV pra ver se encontro essa mesma diferença. De todas as formas, obrigado pela ajuda! :) Commented 14/12/2013 às 16:29
  • Disponha! Eu entendo pouco de matemática, mas também me ocorreu a possibilidade de que sua função tenha uma descontinuidade justamente nesse ponto, o que faria com que o resultado ficasse radicalmente diferente para o valor de 3 ou para um valor tendendo a 3. Não sei se se aplica ao seu caso em particular ou não, mas se o seu teste com C++ der resultados semelhantes vale a pena investigar essa possibilidade...
    – mgibsonbr
    Commented 14/12/2013 às 22:59
  • Entendo, mas não sei se é o caso porque trata-se de uma função senoidal em que o 3 é o comprimento da onda (wavelength, ou lambda). Segundo várias fontes que eu li, nesse domínio (processamento de imagens) o parâmetro lambda é dado em pixels. Por isso, eu não esperava ver muita diferença das ondas produzidas pela imagem 1. Commented 14/12/2013 às 23:08
  • 1
    @mgibsonbr: Fiz o teste em C++ com OpenCV. Lá, os resultados são consistentes, isto é, não faz absolutamente nenhuma diferença usar 3.0 ou 3.00001 no valor de lambda. Creio que a diferença deve ser algum detalhe do pacote scikit-image mesmo, talvez algo como o que você sugeriu. Commented 18/12/2013 às 17:23

2 Respostas 2

9
+50

As diferenças decorreram na verdade de uma pequena confusão (de minha parte mesmo, desculpem!) a respeito do primeiro parâmetro utilizado na chamada da função gabor_kernel do pacote scikit-image. NÃO HÁ ERRO no código de gabor_kernel, o scikit-image apenas utiliza uma abordagem diferente do OpenCV (com o qual eu estava acostumado).

Essa resposta corrige a anterior, e foi alcançada com a ajuda dos colegas em comentários efetuados na questão e também com uma dúvida enviada diretamente ao desenvolvedor do pacote scikit-image (veja issue no Git do projeto neste link).

Explicando melhor...

O filtro de Gabor é composto de um sinal sinusoidal e um envelope gaussiano. O sinal sinuisoidal é definido por três parâmetros: um relacionado ao comprimento da onda (ou à sua frequência, dependendo do ponto de vista preferido), outro à orientação espacial do sinal (rotação) e o terceiro ao deslocamento (offset) da função. O envelope gaussiano é definido por dois parâmetros: um desvio padrão (a média é o centro do kernel) e uma razão de aspecto (quão circular a gaussiana bidimensional é; 1 significa totalmente circular - o scikit-image utiliza dois desvios padrões distintos para x e y, para ter o mesmo efeito).

Talvez seja mais comum na literatura (vide equações disponíveis na Wikipedia, implementação do OpenCV e applet Java referida na resposta anterior) a caracterização do sinal sinusoidal em termos do comprimento de onda (lambda). De todas as formas, o comprimento e a frequência da onda se relacionam segundo uma razão (vide esse link na Wikipedia, por exemplo) tal que a frequencia f é igual a uma velocidade v dividida pelo comprimento de onda lambda:

inserir a descrição da imagem aqui

Assumindo-se 1 para a velocidade (devido ao contexto de processamento de imagens digitais), é uma alternativa totalmente correta calcular o componente complexo do filtro de Gabor com a multiplicação do valor da FREQUÊNCIA ao invés da divisão do valor do COMPRIMENTO DE ONDA.

A divisão que eu mencionei como "DÚVIDA 1" na resposta anterior deve-se apenas à normalização do envelope gaussiano (para que a soma dos pesos no kernel gaussiano seja igual a 1 - notem nos testes abaixo que a escala de valores no eixo z é menor do que nos testes que eu fiz na resposta anterior).

Concluindo...

Na minha construção eu estava considerando que o parâmetro inicial da função gabor_kernel era o comprimento de onda em pixels, e não a frequência do sinal como esperado pelo pacote scikit-image. Desse modo, as variações observadas não eram decorrentes de diferenças de precisão no Python ou de erros no pacote scikit-image, mas sim do fato de que há diferenças significativas entre 3 e 3.00001 em termos de representação de frequência.

Considerando isso, eu mantive no meu código original a definição do kernel em termos do comprimento de onda (lambda), e apenas utilizei o valor 1.0 / fLambda como parâmetro para a chamada de gabor_kernel (isto é, simplesmente converto o valor do comprimento de onda para o valor de frequência, como esperado pelo pacote).

Testando novamente, o kernel agora apresenta os resultados consistentes e sem variações. :)

Teste com fLambda = 3: inserir a descrição da imagem aqui

Teste com fLambda = 3.00001: inserir a descrição da imagem aqui

Código final correto:

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np
from scipy import interpolate
from mpl_toolkits.mplot3d import Axes3D
from pylab import *
from skimage.filter import gabor_kernel
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

# Função utilizada para reescalar o kernel
def resize_kernel(aKernelIn, iNewSize):
    x = np.array([v for v in range(len(aKernelIn))])
    y = np.array([v for v in range(len(aKernelIn))])
    z = aKernelIn

    xx = np.linspace(x.min(), x.max(), iNewSize)
    yy = np.linspace(y.min(), y.max(), iNewSize)

    aKernelOut = np.zeros((iNewSize, iNewSize), np.float)    
    oNewKernel = interpolate.RectBivariateSpline(x, y, z)
    aKernelOut = oNewKernel(xx, yy)

    return aKernelOut


if __name__ == "__main__":
    fLambda = 3              # comprimento de onda (em pixels)
    fTheta  = 0              # orientação (em radianos)
    fSigma  = 0.56 * fLambda # envelope gaussiano (com 1 oitava de largura de banda)
    fPsi    = np.pi / 2      # deslocamento (offset)

    # Tamanho do kernel (3 desvios para cada lado, para limitar cut-off)
    iLen = int(math.ceil(3.0 * fSigma))
    if(iLen % 2 != 0):
        iLen += 1

    # Obtém o kernel Gabor para os parâmetros definidos
    z = np.real(gabor_kernel(1.0/fLambda, theta=fTheta, sigma_x=fSigma, sigma_y=fSigma, offset=fPsi))

    # Plotagem do kernel

    fig = plt.figure(figsize=(16, 9))
    fig.suptitle(r'Gabor kernel para $\lambda=3.0$, $\theta=0.0$, $\sigma=0.56\lambda$ e $\psi=\frac{\pi}{2}$', fontsize=25)
    grid = gridspec.GridSpec(1, 2, width_ratios=[1, 2])

    # Gráfico 2D
    plt.gray()
    ax = fig.add_subplot(grid[0])
    ax.set_title(u'Visualização 2D')
    ax.imshow(z, interpolation='bicubic')
    ax.set_xticklabels([v for v in range((-iLen/2)-1, iLen/2+1)], fontsize=10)
    ax.set_yticklabels([v for v in range((-iLen/2)-1, iLen/2+1)], fontsize=10)

    # Gráfico em 3D

    # Reescalona o kernel para uma exibição melhor
    z = resize_kernel(z, 300)
    # Eixos x e y no intervalo do tamanho do kernel
    x = np.linspace(-iLen/2, iLen/2, 300)
    y = x
    x, y = meshgrid(x, y)

    ax = fig.add_subplot(grid[1], projection='3d')
    ax.set_title(u'Visualização 3D')
    ax.plot_surface(x, y, z, cmap='hot')

    print(iLen)

    plt.show()
7

Com a ajuda dos comentários efetuados, me parece que há um possível erro no cálculo realizado pela função gabor_kernel do pacote scikit-image.

A função original (obrigado ao @ricidleiv por postar o link com o código) faz o cálculo do valor complexo. As linhas que executam tal cálculo estão no final da função, e são reproduzidas abaixo:

g = np.zeros(y.shape, dtype=np.complex)
g[:] = np.exp(-0.5 * (rotx**2 / sigma_x**2 + roty**2 / sigma_y**2))
g /= 2 * np.pi * sigma_x * sigma_y                                   # <- DÚVIDA 1
g *= np.exp(1j * (2 * np.pi * frequency * rotx + offset))            # <- DÚVIDA 2

return g

Para comparação, considere o código da mesma função na implementação do OpenCV (disponível aqui) e as fórmulas (em versão complexa, real e imaginária) disponíveis na Wikipedia:

inserir a descrição da imagem aqui

A implementação do OpenCV calcula diretamente o valor real (que é o que me interessa, e de fato a implementação que eu fiz no Python chama np.real para extrair o valor real do valor complexo devolvido). A implementação do scikit-image calcula o valor complexo, mas a implementação tem algumas características que parecem estar incorretas.

Na linha marcada com o comentário "<- DÚVIDA 1", por exemplo, por que essa divisão é executada? Ela parece ser um valor redundante, uma vez que a divisão por 2 * sigma da fórmula já está incluida pela multiplicação do valor -0.5 na linha anterior.

Além disso, a parte final da fórmula - que calcula o sinal complexo, está multiplicando a frequencia (lambda, na fórmula) pelo valor de x rotacionado, sendo que deveria dividi-lo ao meu entender.

Eu fiz um pequeno teste copiando a função original do scikit-image para o meu código e chamando-a my_gabor_kernel. Alterei essa parte do código de forma que me parecesse correto segundo a fórmula do filtro Gabor na Wikipedia. Ficou assim:

g = np.zeros(y.shape, dtype=np.complex)
g[:] = np.exp(-0.5 * (rotx**2 / sigma_x**2 + roty**2 / sigma_y**2))
g *= np.exp(1j * (2 * np.pi * rotx / frequency + offset))

return g

Após executar novos testes, não há mais diferença entre a utilização do valor 3 ou 3.0001. O resultado do kernel produzido também parece condizente, ainda que com uma frequencia aparente (visual) leventemente maior do que anteriormente observado. Considerando que o valor de lambda é dado em pixels, esse resultado me parece condizente e o "problema" de precisão não ocorre.

Usando o valor 3: Usando o valor 3

Usando o valor 3.000001: Usando o valor 3.0000001

EDIT:

Aqui está o código completo que eu utilizei no meu teste, caso alguém mais deseje testar:

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np
from scipy import interpolate
from mpl_toolkits.mplot3d import Axes3D
from pylab import *
from skimage.filter import gabor_kernel
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

def my_gabor_kernel(frequency, theta=0, bandwidth=1, sigma_x=None, sigma_y=None,
                 offset=0):
    if sigma_x is None:
        sigma_x = _sigma_prefactor(bandwidth) / frequency
    if sigma_y is None:
        sigma_y = _sigma_prefactor(bandwidth) / frequency

    n_stds = 3
    x0 = np.ceil(max(np.abs(n_stds * sigma_x * np.cos(theta)),
                     np.abs(n_stds * sigma_y * np.sin(theta)), 1))
    y0 = np.ceil(max(np.abs(n_stds * sigma_y * np.cos(theta)),
                     np.abs(n_stds * sigma_x * np.sin(theta)), 1))
    y, x = np.mgrid[-y0:y0+1, -x0:x0+1]

    rotx = x * np.cos(theta) + y * np.sin(theta)
    roty = -x * np.sin(theta) + y * np.cos(theta)

    g = np.zeros(y.shape, dtype=np.complex)
    g[:] = np.exp(-0.5 * (rotx**2 / sigma_x**2 + roty**2 / sigma_y**2))
    #g /= 2 * np.pi * sigma_x * sigma_y
    g *= np.exp(1j * (2 * np.pi * rotx / frequency + offset))

    return g

# Função utilizada para reescalar o kernel
def resize_kernel(aKernelIn, iNewSize):
    x = np.array([v for v in range(len(aKernelIn))])
    y = np.array([v for v in range(len(aKernelIn))])
    z = aKernelIn

    xx = np.linspace(x.min(), x.max(), iNewSize)
    yy = np.linspace(y.min(), y.max(), iNewSize)

    aKernelOut = np.zeros((iNewSize, iNewSize), np.float)    
    oNewKernel = interpolate.RectBivariateSpline(x, y, z)
    aKernelOut = oNewKernel(xx, yy)

    return aKernelOut


if __name__ == "__main__":
    fLambda = 3.0000001      # comprimento de onda (em pixels)
    fTheta  = 0              # orientação (em radianos)
    fSigma  = 0.56 * fLambda # envelope gaussiano (com 1 oitava de largura de banda)
    fPsi    = np.pi / 2      # deslocamento (offset)

    # Tamanho do kernel (3 desvios para cada lado, para limitar cut-off)
    iLen = int(math.ceil(3.0 * fSigma))
    if(iLen % 2 != 0):
        iLen += 1

    # Obtém o kernel Gabor para os parâmetros definidos
    z = np.real(my_gabor_kernel(fLambda, theta=fTheta, sigma_x=fSigma, sigma_y=fSigma, offset=fPsi))
    print(z.shape)

    # Plotagem do kernel

    fig = plt.figure(figsize=(16, 9))
    fig.suptitle(r'Gabor kernel para $\lambda=3.0$, $\theta=0.0$, $\sigma=0.56\lambda$ e $\psi=\frac{\pi}{2}$', fontsize=25)
    grid = gridspec.GridSpec(1, 2, width_ratios=[1, 2])

    # Gráfico 2D
    plt.gray()
    ax = fig.add_subplot(grid[0])
    ax.set_title(u'Visualização 2D')
    ax.imshow(z, interpolation='bicubic')
    ax.set_xticklabels([v for v in range((-iLen/2)-1, iLen/2+1)], fontsize=10)
    ax.set_yticklabels([v for v in range((-iLen/2)-1, iLen/2+1)], fontsize=10)

    # Gráfico em 3D

    # Reescalona o kernel para uma exibição melhor
    z = resize_kernel(z, 300)
    # Eixos x e y no intervalo do tamanho do kernel
    x = np.linspace(-iLen/2, iLen/2, 300)
    y = x
    x, y = meshgrid(x, y)

    ax = fig.add_subplot(grid[1], projection='3d')
    ax.set_title(u'Visualização 3D')
    ax.plot_surface(x, y, z, cmap='hot')

    print(iLen)

    plt.show()

P.S.: Também para comparação, pode-se brincar com essa applet Java do depto de computação da Universidade de Groningen. Pra chegar em um resultado semelhante, use 90 no deslocamento (a applet recebe o offset em graus). Como a imagem tem tamanho fixo (o meu código escala pra facilitar a visualização), use 30 ao invés de 3 no comprimento para ser mais fácil de enxergar. Exemplo do uso da applet:

inserir a descrição da imagem aqui

EDIT:

Eu abri uma issue no Git do projeto scikit-image (link aqui). Quando tiver uma resposta mais oficial eu atualizo aqui.

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .