Graph differences in matplotlib for different parameter accuracies

The differences were actually a little confusion (for my part, sorry!) regarding the first parameter used in the call of the gabor_kernel function of the scikit-image package. THERE IS NO ERROR in the gabor_kernel code, scikit-image only uses an approach other than OpenCV (with which I was accustomed).

This response corrects the above, and was reached with the help of colleagues in comments made on the issue and also with a question sent directly to the developer of the scikit-image package (see issue in the project's Git in this link ).

Explaining better ...

The Gabor filter is composed of a sinusoidal signal and a Gaussian envelope. The sinusoidal signal is defined by three parameters: one related to the wavelength (or its frequency depending on the preferred point of view), the other to the spatial orientation of the signal (rotation) and the third to the offset of the function. The Gaussian envelope is defined by two parameters: a standard deviation (the mean is the center of the kernel) and an aspect ratio (how circular the two-dimensional Gaussian is; 1 means totally circular - the scikit-image uses two distinct standard deviations for x and y , to have the same effect).

Perhaps it is more common in the literature (see equations available in the Wikipedia , OpenCV implementation and Java applet referenced in the previous answer ) the characterization of the sinusoidal signal in terms of the wavelength (lambda). However, the length and frequency of the wave are related to one reason (see this link on Wikipedia , for example) such that the frequency f is equal to a velocity v divided by the lambda wavelength:

Assuming1forspeed(duetothecontextofdigitalimageprocessing),itisaentirelycorrectalternativetocalculatethecomplexcomponentoftheGaborfilterbymultiplyingtheFREQUENCYvalueinsteadofthedivisionofthevalueoftheWAVELENGTH.

ThedivisionImentionedas"DOUBLE 1" in the previous answer is due only to the normalization of the Gaussian envelope (so that the sum of the weights in the Gaussian kernel equals 1 - note in the tests below that the scale of values in z axis is less than in the tests I did in the previous answer).

Finishing ...

In my construction I was considering that the initial parameter of the gabor_kernel function was the wavelength in pixels, not the signal frequency as expected by the scikit-image package. Thus, the observed variations were not due to differences in Python precision or errors in the scikit-image packet, but to the fact that there are significant differences between 3 and 3,00001 in terms of frequency representation.

Considering this, I kept the kernel definition in terms of the wavelength (lambda) in my original code, and just used the value 1.0 / fLambda as the parameter for the gabor_kernel call is, simply convert the value of the wavelength to the frequency value, as expected by the packet).

Testing again, the kernel now presents consistent, unvarying results. :)

Test with fLambda = 3:

TestwithfLambda=3.00001:

Correct end code:

#!/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()

With the help of the comments made, it seems to me that there is a possible error in the calculation performed by the gabor_kernel function of the scikit-image package.

The original function (thanks to @ricidleiv for posting the link with the code) computes the complex value. Lines that perform such a calculation are at the end of the function, and are reproduced below:

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

For comparison, consider the code for the same function in the OpenCV implementation (available here ) and the formulas (in complex, real and imaginary version) available at Wikipedia :

TheOpenCVimplementationdirectlycalculatestheactualvalue(whichiswhatinterestsme,andinfacttheimplementationIdidinPythoncallsnp.realtoextracttheactualvalueofthereturnedcomplexvalue).Thescikit-imageimplementationcalculatesthecomplexvalue,buttheimplementationhassomefeaturesthatappeartobeincorrect.

Onthelinemarkedwiththecomment"< - DOUBT 1", for example, why is this division performed? It appears to be a redundant value since the 2 * sigma division of the formula is already included by multiplying the value -0.5 on the previous line.

In addition, the final part of the formula - which computes the complex signal, is multiplying the frequency (lambda, in the formula) by the value of x rotated, and should divide it as I understand it.

I did a little test copying the original scikit-image function to my code and calling it my_gabor_kernel. I changed that part of the code so that it would look right to me according to the Gabor filter formula on Wikipedia. It looks like this:

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

After running new tests, there is no difference between using the 3 or 3.0001 value. The result of the kernel produced also seems to be consistent, albeit with an apparent (visual) frequency slightly higher than previously observed. Considering that the lambda value is given in pixels, this result seems to me to be correct and the precision "problem" does not occur.

Using the value 3:

Usingthe3.000001value:

EDIT:

Here is the complete code I used in my test if someone else wants to test:

#!/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()

PS: Also for comparison, one can play with this Java applet from the computing depot of the University of Groningen. To arrive at a similar result, use 90 in the offset (the applet receives the offset in degrees). Since the image has a fixed size (my code scales for easy viewing), use 30 instead of 3 in length to make it easier to see. Example of using the applet:

EDIT:

I'veopenedaissueinGit'sscikit-imageproject( link here ). When I have a more official answer I update here.