Plot differences in matplotlib for different parameter accuracies

Question

For a course work, I built the following code in Python using the Matplotlib and scikit-image packages:

#!/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 the definition of the variable fLambda in line 29:

fLambda = 3.0000000001

If the variable is defined in this way (i.e., with a large number of decimals, but with a value close to 0), the Gabor kernel plotting (which is the intention of this program) appears as expected:

However, if the variable is defined in one of the forms below (i.e., with the round value '3'):

fLambda = 3.0
# ou
fLambda = float(3)

[...] plotting appears with a quite different result than expected:

My first impression was that the interpolation used in the kernel rescaling (resize_kernel function) was wrong, but can be seen in the 2D view (which directly uses the kernel returned by the scikit-image package, i.e., NO interpolation)that the kernel is already different.

Other information: when using the value 3.000000000000001 (15 decimals) the result is equal to the first image; when using values with more decimals (for example, 3.000000000001 - 1 decimal place), the result is already equal to the second image.

Does this result seem right to anyone (is the Gabor kernel so different for integer and floating point values)? Might there be some accuracy error involved? Is it some Python specificity or even scikit-image package?

Python version: 2.7.5 (windows 7 32bit) Version of Matplotlib: 1.3.1 Version of numpy: 1.8.0 Scikit-image version: 0.9.3

Answer 1

The differences were in fact a little confusion (on my own, sorry! ) regarding the first parameter used in the gabor_kernel function call of the scikit-image package. THERE IS NO ERROR in gabor_kernel code, scikit-image only uses a different approach to Opencv (which I was used to).

This answer corrects the previous one, and was achieved with the help of colleagues in comments made in the question and also with a question sent directly to the scikit-image package developer (see Git Content of the project at this link).

Explaining better...

The Gabor filter is composed of a sinusoidal signal and a Gaussian envelope. The synuisoidal 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 displacement (offset) of the function. The Gaussian envelope is defined by two parameters: a standard deviation (the mean is the kernel center) and an aspect ratio (how circular the two-dimensional Gaussian is; 1 means fully 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 Wikipedia, implementation of the Opencv and Java applet referred to in the previous response) characterization of the sinusoidal signal in terms of wavelength (lambda). In all forms, the length and frequency of the wave are related according to a ratio (see this link on Wikipedia, for example) such that the frequency f is equal to a speed v divided by the wavelength lambda:

Assuming 1 for speed (due to the context of digital image processing), it is an alternative fully correct calculate the complex component of the Gabor filter by multiplying the FREQUENCY value instead of dividing the WAVELENGTH value.

The division I mentioned as "DOUBT 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 is equal to 1 - note in the tests below that the scale of values in the z axis is smaller than in the tests I did in the previous answer).

Concluding...

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

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

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

Test with fLambda = 3:

Test with fLambda = 3.00001:

Correct final 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()

Answer 2

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 link with the code) makes the calculation of the complex value. The lines that perform such 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 code from the same function in the Opencv implementation (available here) and the formulas (in complex, real and imaginary versions) available in Wikipedia:

The Opencv implementation directly calculates the real value (which is what interests me, and in fact the implementation I made in Python calls np.real to extract the real value from the returned complex value). The implementation of scikit-image calculates the complex value, but the implementation has some features that seem to be incorrect.

On the line marked with the comment "<- DOUBT 1", for example, why is this division executed? It seems to be a redundant value, since the division by 2 * sigma of the formula is already included by multiplying the value -0.5 in the previous line.

Also, the final part of the formula - which calculates the complex signal, is multiplying the frequency (lambda, in the formula) by the value of x rotated, and should divide it to my mind.

I did a little test by copying the original scikit-image function to my code and calling it my_gabor_kernel. I changed this part of the code so that it seemed correct according to the formula of the Gabor filter on Wikipedia. It was 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 longer any difference between using the value 3 or 3.0001. The result of the produced kernel also seems 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 suitable and the precision "problem" does not occur.

Using the value 3:

Using the value 3.000001:

EDIT:

Here’s the full code I used in my test, in case anyone 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()

P.S.: Also for comparison, one can play with this Java applet from the Computer Department at 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 scale code for easy viewing), use 30 instead of 3 in length to make it easier to see. Example of using the applet:

EDIT:

I opened an Issue in Git from the scikit-image project (link here). When you have a more official reply I will update here.