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:
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:
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
3.0000000001
é maior que3
, se em qualquer momento for exigido um teto ele será4
no primeiro caso e3
no segundo (quando vi oceil
no seu código quase achei que havia achado o erro, mas aí percebi que ele envolviafSigma
e nãofLambda
). 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 ser3
(e, portanto, o resultado ser igual ao da segunda imagem).3
ou para um valor tendendo a3
. 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...