Gabor Filter implementation in Frequency domain

Question

Here we have the Spatial domain implementation of Gabor filter. But, I need to implement a Gabor filter in the Frequency Domain for performance reasons.

I have found the Frequency Domain equation of Gabor Filter:

^{I am actually in doubt about the correctness and/or applicability of this formula.}

Source Code

So, I have implemented the following :

public partial class GaborFfftForm : Form
{
    private double Gabor(double u, double v, double f0, double theta, double a, double b)
    {
        double rad = Math.PI / 180 * theta;
        double uDash = u * Math.Cos(rad) + v * Math.Sin(rad);
        double vDash = (-1) * u * Math.Sin(rad) + v * Math.Cos(rad);

        return Math.Exp((-1) * Math.PI * Math.PI * ((uDash - f0) / (a * a)) + (vDash / (b * b)));
    }

    public Array2d<Complex> GaborKernelFft(int sizeX, int sizeY, double f0, double theta, double a, double b)
    {
        int halfX = sizeX / 2;
        int halfY = sizeY / 2;

        Array2d<Complex> kernel = new Array2d<Complex>(sizeX, sizeY);

        for (int u = -halfX; u < halfX; u++)
        {
            for (int v = -halfY; v < halfY; v++)
            {
                double g = Gabor(u, v, f0, theta, a, b);

                kernel[u + halfX, v + halfY] = new Complex(g, 0);
            }
        }

        return kernel;
    }

    public GaborFfftForm()
    {
        InitializeComponent();

        Bitmap image = DataConverter2d.ReadGray(StandardImage.LenaGray);
        Array2d<double> dImage = DataConverter2d.ToDouble(image);

        int newWidth = Tools.ToNextPowerOfTwo(dImage.Width) * 2;
        int newHeight = Tools.ToNextPowerOfTwo(dImage.Height) * 2;

        double u0 = 0.2;
        double v0 = 0.2;
        double alpha = 10;//1.5;
        double beta = alpha;

        Array2d<Complex> kernel2d = GaborKernelFft(newWidth, newHeight, u0, v0, alpha, beta);

        dImage.PadTo(newWidth, newHeight);
        Array2d<Complex> cImage = DataConverter2d.ToComplex(dImage);
        Array2d<Complex> fImage = FourierTransform.ForwardFft(cImage);

        // FFT convolution .................................................
        Array2d<Complex> fOutput = new Array2d<Complex>(newWidth, newHeight);
        for (int x = 0; x < newWidth; x++)
        {
            for (int y = 0; y < newHeight; y++)
            {
                fOutput[x, y] = fImage[x, y] * kernel2d[x, y];
            }
        }

        Array2d<Complex> cOutput = FourierTransform.InverseFft(fOutput);
        Array2d<double> dOutput = Rescale2d.Rescale(DataConverter2d.ToDouble(cOutput));

        //dOutput.CropBy((newWidth-image.Width)/2, (newHeight - image.Height)/2);

        Bitmap output = DataConverter2d.ToBitmap(dOutput, image.PixelFormat);

        Array2d<Complex> cKernel = FourierTransform.InverseFft(kernel2d);
        cKernel = FourierTransform.RemoveFFTShift(cKernel);
        Array2d<double> dKernel = DataConverter2d.ToDouble(cKernel);
        Array2d<double> dRescaledKernel = Rescale2d.Rescale(dKernel);
        Bitmap kernel = DataConverter2d.ToBitmap(dRescaledKernel, image.PixelFormat);

        pictureBox1.Image = image;
        pictureBox2.Image = kernel;
        pictureBox3.Image = output;
    }
}

^{Just concentrate on the algorithmic steps at this time.}

I have generated a Gabor kernel in the frequency domain. Since, the kernel is already in Frequency domain, I didn't apply FFT to it, whereas image is FFT-ed. Then, I multiplied the kernel and the image to achieve FFT-Convolution. Then they are inverse-FFTed and converted back to Bitmap as usual.

Output

1. The kernel looks okay. But, The filter-output doesn't look very promising (or, does it?).
2. The orientation (theta) doesn't have any effect on the kernel.
3. The calculation/formula is frequently suffering from divide-by-zero exception up on changing values.

How can I fix those problems?

Oh, and, also,

• what do the parameters α, β, represent?
• what should be the appropriate value of f₀?

---

Update:

I have modified my code according to @Cris Luoengo's answer.

    private double Gabor(double u, double v, double u0, double v0, double a, double b)
    {
        double p = (-2) * Math.PI * Math.PI;
        double q = (u-u0)/(a*a);
        double r = (v - v0) / (b * b);

        return Math.Exp(p * (q + r));
    }

    public Array2d<Complex> GaborKernelFft(int sizeX, int sizeY, double u0, double v0, double a, double b)
    {
        double xx = sizeX;
        double yy = sizeY;
        double halfX = (xx - 1) / xx;
        double halfY = (yy - 1) / yy;

        Array2d<Complex> kernel = new Array2d<Complex>(sizeX, sizeY);

        for (double u = 0; u <= halfX; u += 0.1)
        {
            for (double v = 0; v <= halfY; v += 0.1)
            {
                double g = Gabor(u, v, u0, v0, a, b);

                int x = (int)(u * 10);
                int y = (int)(v * 10);

                kernel[x,y] = new Complex(g, 0);
            }
        }

        return kernel;
    }

where,

        double u0 = 0.2;
        double v0 = 0.2;
        double alpha = 10;//1.5;
        double beta = alpha;

I am not sure whether this is a good output.

Answer 1

There seems to be a typo in the equation for the Gabor filter that you found. The Gabor filter is a translated Gaussian in the frequency domain. Hence, it needs to have u² and v² in the exponent.

Equation (2) in your link seems more sensible, but still misses a 2:

exp( -2(πσ)² (u-f₀)² )

It is the 1D case, this is the filter we want to use in the direction θ. We now multiply in the perpendicular direction, v, with a non-shifted Gaussian. I set α and β to be the inverse of the two sigmas:

exp( -2(π/α)² (u-f₀)² ) exp( -2(π/β)² v² ) = exp( -2π²((u-f₀)/α)² + -2π²(v/β)² )

You should implement the above equation with u and v rotated over θ, as you already do.

Also, u and v should run from -0.5 to 0.5, not from -sizeX/2 to sizeX/2. And that is assuming your FFT sets the origin in the middle of the image, which is not common. Typically the FFT algorithms set the origin in a corner of the image. So you should probably have your u and v run from 0 to (sizeX-1)/sizeX instead.

With a corrected implementation as above, you should set f₀ to be between 0 and 0.5 (try 0.2 to start with), and α and β should be small enough such that the Gaussian doesn't reach the 0 frequency (you want the filter to be 0 there)

In the frequency domain, your filter will look like a rotated Gaussian away from the origin.

In the spatial domain, the magnitude of your filter should look again like a Gaussian. The imaginary component should look like this (picture links to Wikipedia page I found it on):

(i.e. it is anti-symmetric (odd) in the θ direction), possibly with more lobes depending on α, β and f₀. The real component should be similar but symmetric (even), with a maximum in the middle. Note that after IFFT, you might need to shift the origin from the top-left corner to the middle of the image (Google "fftshift").

---

Note that if you set α and β to be equal, the rotation of the u-v plane is irrelevant. In this case, you can use cartesian coordinates instead of polar coordinates to define the frequency. That is, instead of defining f₀ and θ as parameters, you define u₀ and v₀. In the exponent you then replace u-f₀ with u-u₀, and v with v-v₀.

---

The code after the edit of the question misses the square again. I would write the code as follows:

private double Gabor(double u, double v, double u0, double v0, double a, double b)
{
    double p = (-2) * Math.PI * Math.PI;
    double q = (u-u0)/a;
    double r = (v - v0)/b;

    return Math.Exp(p * (q*q + r*r));
}

public Array2d<Complex> GaborKernelFft(int sizeX, int sizeY, double u0, double v0, double a, double b)
{
    double halfX = sizeX / 2;
    double halfY = sizeY / 2;

    Array2d<Complex> kernel = new Array2d<Complex>(sizeX, sizeY);

    for (double y = 0; y < sizeY; y++)
    {
        double v = y / sizeY;
        // v -= HalfY;  // whether this is necessary or not depends on your FFT
        for (double x = 0; x < sizeX; x++)
        {
            double u = x / sizeX;
            // u -= HalfX;  // whether this is necessary or not depends on your FFT
            double g = Gabor(u, v, u0, v0, a, b);

            kernel[(int)x, (int)y] = new Complex(g, 0);
        }
    }

    return kernel;
}