Gabor filter implementation

Question

I have implemented the Gabor Filter of which I provided the source code as bellow.

The code is deadly slow.

Applying the filter on a grayscale image of Lena 512x512 takes 7577 milliseconds.

And the filter bank (12 filters) freezes and doesn't even display any value.

Source Code:

Here, you can obtain the complete source code as a zip file.

Here are the most relevant parts of the source code,

Gabor kernel,

public class GaborKernel
{
    private Matrix __gaborKernelData;
    public GaborKernelKind GaborKernelKind { get; set; }
    public int Width { get { return __gaborKernelData.Rows; } }
    public int Height { get { return __gaborKernelData.Cols; } }
    public double Size { get; set; }
    public double Lambda { get; set; }
    public double Theta { get; set; }
    public double Psi { get; set; }
    public double Sigma { get; set; }
    public double Gamma { get; set; }
    public bool Normalized { get; set; }

    public GaborKernel() {}

    public GaborKernel(double size,double lambda,double theta,double psi,double sigma,double gamma, bool normalized, GaborKernelKind gaborKernelType)
    {
        Size = size;
        Lambda = lambda;
        Theta = theta;
        Psi = psi;
        Sigma = sigma;
        Gamma = gamma;
        Normalized = normalized;
        GaborKernelKind = gaborKernelType;
    }

    public double this[int x, int y]
    {
        get
        {
            return __gaborKernelData[x,y];
        }
        set
        {
            __gaborKernelData[x,y] = value;
        }
    }

    public void Compute()
    {
        double sigmaX = Sigma;
        double sigmaY = Sigma / Gamma;

        double a = Math.Max(Math.Abs(Size * sigmaX * Math.Cos(Theta)), Math.Abs(Size * sigmaY * Math.Sin(Theta)));
        int xMax = (int) Math.Ceiling(Math.Max(1, a));

        double b = Math.Max(Math.Abs(Size * sigmaX * Math.Sin(Theta)), Math.Abs(Size * sigmaY * Math.Cos(Theta)));
        int yMax = (int)Math.Ceiling(Math.Max(1, b));

        List<int> xValues = Matrix.Vector(-xMax, xMax, increment: 1);
        List<int> yValues = Matrix.Vector(-yMax, yMax, increment: 1);

        System.Diagnostics.Debug.Assert(xValues.Count == (2 * xMax + 1));
        System.Diagnostics.Debug.Assert(yValues.Count == (2 * yMax + 1));

        __gaborKernelData = new Matrix(xValues.Count, yValues.Count);

        double sum = 0;

        switch (GaborKernelKind)
        {
            case GaborKernelKind.Real:
                for (int i = 0; i < xValues.Count; i++)
                {
                    for (int j = 0; j < yValues.Count; j++)
                    {
                        sum += __gaborKernelData[i, j] = GaborFunction.RealFunction2D(xValues[i], yValues[j], Lambda, Theta, Psi, Sigma, Gamma);
                    }
                }
                break;

            case GaborKernelKind.Imaginary:
                for (int i = 0; i < xValues.Count; i++)
                {
                    for (int j = 0; j < yValues.Count; j++)
                    {
                        sum += __gaborKernelData[i, j] = GaborFunction.ImaginaryFunction2D(xValues[i], yValues[j], Lambda, Theta, Psi, Sigma, Gamma);
                    }
                }
                break;

            case GaborKernelKind.Magnitude:
                for (int i = 0; i < xValues.Count; i++)
                    for (int j = 0; j < yValues.Count; j++)
                        sum += __gaborKernelData[i, j] = GaborFunction.Function2D(xValues[i], yValues[j], Lambda, Theta, Psi, Sigma, Gamma).Magnitude;
                string str3 = string.Empty;
                break;

            case GaborKernelKind.SquaredMagnitude:
                for (int i = 0; i < xValues.Count; i++)
                    for (int j = 0; j < yValues.Count; j++)
                    {
                        double v = GaborFunction.Function2D(xValues[i], yValues[j], Lambda, Theta, Psi, Sigma, Gamma).Magnitude;
                        sum += __gaborKernelData[i, j] = v * v;
                    }
                string str4 = string.Empty;
                break;
        }

        if (Normalized)
        {
            __gaborKernelData = __gaborKernelData / sum;
        }
    }
}

Gabor filter,

public class GaborFilter
{
    public GaborKernel GaborKernel { get; set; }
    public int GaborKernelSize { get; set; }
    public double Theta { get; set; }
    public double Lambda { get; set; }
    public double Psi { get; set; }
    public double Sigma { get; set; }
    public double Gamma { get; set; }
    public double F
    {
        get { return 1 / Lambda; }
        set { Lambda = 1/value;} 
    }

    public void Process()
    {
        GaborKernel = new GaborKernel(GaborKernelSize, Lambda, Theta, Psi, Sigma, Gamma, true, GaborKernelKind.Imaginary);
        GaborKernel.Compute();
    }

    public GaborFilter(){}

    public Bitmap Apply(Bitmap srcImage)
    {
        if (GaborKernel != null)
        {
            // lock source bitmap data
            BitmapData srcData = srcImage.LockBits(
                new Rectangle(0, 0, srcImage.Width, srcImage.Height),
                ImageLockMode.ReadOnly, srcImage.PixelFormat);

            Bitmap dstImage = Grayscale.CreateGrayscaleImage(srcImage.Width, srcImage.Height);

            BitmapData dstData = dstImage.LockBits(new Rectangle(0, 0, srcImage.Width, srcImage.Height),
                                    ImageLockMode.ReadOnly, PixelFormat.Format8bppIndexed);

            try
            {
                // apply the filter
                Apply(srcData, dstData);

                if ((srcImage.HorizontalResolution > 0) && (srcImage.VerticalResolution > 0))
                {
                    dstImage.SetResolution(srcImage.HorizontalResolution, srcImage.VerticalResolution);
                }
            }
            finally
            {
                // unlock source image
                srcImage.UnlockBits(srcData);
                dstImage.UnlockBits(dstData);
            }

            return dstImage;
        }
        else
        {
            throw new Exception("Kernel is not processed yet!");
        }
    }

    unsafe private void Apply(BitmapData sourceData, BitmapData destinationData)
    {
        int kernelHeight = GaborKernel.Width;
        int kernelWidth = GaborKernel.Height;

        int centerX = kernelHeight / 2;
        int centerY = kernelWidth / 2;

        int width = sourceData.Width;
        int height = sourceData.Height;

        int srcStride = sourceData.Stride;
        int srcOffset = srcStride - width;

        byte* address = (byte*)sourceData.Scan0.ToPointer();


        int[,] response = new int[height, width];

        int max = int.MinValue;
        int min = int.MaxValue;

        // for each image row
        for (int y = 0; y < height; y++)
        {
            // for each pixel in the row
            for (int x = 0; x < width; x++, address++)
            {
                double sum = 0;

                // for each kernel row
                for (int i = 0; i < kernelHeight; i++)
                {
                    int ir = i - centerY;
                    int t = y + ir;

                    // skip row
                    if (t < 0)
                        continue;

                    // break
                    if (t >= height)
                        break;

                    int col = ir * srcStride;

                    // for each kernel value in the row
                    for (int j = 0; j < kernelWidth; j++)
                    {
                        int jr = j - centerX;
                        t = x + jr;

                        // skip column
                        if (t < 0)
                            continue;

                        if (t < width)
                        {
                            double k = GaborKernel[i, j];
                            sum += k * address[col + jr];
                        }
                    }

                    int v = response[y, x] = (int)sum;

                    if (v > max) max = v;
                    if (v < min) min = v;
                }
            }

            address += srcOffset;
        }

        byte* dst = (byte*)destinationData.Scan0.ToPointer();
        int pixelSize = System.Drawing.Image.GetPixelFormatSize(destinationData.PixelFormat) / 8;
        int dstStride = destinationData.Stride;
        int dstOffset = dstStride - width * pixelSize;

        if (destinationData.PixelFormat == PixelFormat.Format8bppIndexed)
        {
            // for each image row
            for (int y = 0; y < height; y++)
            {
                // for each pixel in the row
                for (int x = 0; x < width; x++, dst++)
                {
                    *dst = (byte)((255 * (response[y, x] - min)) / (max - min));
                }

                dst += dstOffset;
            }
        }
        else
        {
            throw new Exception("Grayscale srcImage needed!");
        }
    }
}

Gabor filter bank,

public class GaborFilterBank
{
    public int NoOfFilters { get; set; }
    public double FilterAngle { get; set; }        
    public int KernelDimension { get; set; }

    public GaborFilterBank() {}

    public List<Bitmap> Apply(Bitmap bitmap)
    {
        List<Bitmap> list = new List<Bitmap>();

        //int size = 3;        // kernel size
        double lambda = 4.0; // wavelength
        //double theta = 0.6;  // orientation
        double psi = 1.0;    // phase offset
        double sigma = 2.0;  // Gaussian variance
        double gamma = 0.3;  // aspect ratio

        ///////////////////////////////////////
        double degrees = FilterAngle;

        GaborFilter filter;

        for (int i = 0; i < NoOfFilters; i++)
        {
            filter = new GaborFilter();
            filter.GaborKernelSize = KernelDimension;
            filter.Lambda = lambda;
            filter.Theta = Tools.ToRadian(FilterAngle);
            filter.Psi = psi;
            filter.Sigma = sigma;
            filter.Gamma = gamma;
            filter.Process();

            Bitmap temp = filter.Apply(bitmap);

            list.Add(temp);

            degrees += FilterAngle;
        }

        return list;
    }
}

Have you tried to run the profiler and see which part is the bottleneck? — t3chb0t, Aug 31 at 16:27
It's a tool in visual studio (or resharper) that can measure a lot of things while your code is running and tell you which method took how long to execute or which one was called most etc. — t3chb0t, Aug 31 at 16:38

Rainer P. · Answer 1 · 2016-08-31 22:07:58Z

Speeding up your code

There's four nested loops in the code and that's obviously the resource hog. As shown in the profiling results, 75% of CPU time is spend in the innermost loop accessing a kernel element via multiple layers of abstraction.

Keep a local copy of the kernel around in form of a simple array, just as you do for the source image. Using unsafe pointers shouldn't be neccessary for either of them.

int[,] img = new int[height, width];
// Fill image
Complex[,] kernel = new Complex[kernelHeight, kernelWidth];
// Fill kernel
Complex[,] response = new Complex[height, width];

for (...

Using complex filters and results should also speed things up a little bit over separate filters for real and imaginary part. It also improves readability and reduces duplicate code, but they're not a must.

In addition, 14% of CPU time is spend in the innermost loops without counting the kernel element access time. Improve it by moving calculations and branches to higher level loops.

for (int imgRow = 0; imgRow < height; imgRow++)
{
    int kernelStartRow = Math.Max(0, 0 - imgRow + centerY);
    int kernelEndRow = Math.Min(kernelHeight, height - imgRow + centerY);
    for (int imgCol = 0; imgCol < width; imgCol++)
    {
        int kernelStartCol = Math.Max(0, 0 - imgCol + centerX);
        int kernelEndCol = Math.Min(kernelWidth, width - imgCol + centerX);    
        Complex sum = 0;
        for (int kernelRow = kernelStartRow; kernelRow < kernelEndRow; kernelRow++)
        {
            for (int kernelCol = kernelStartCol; kernelCol < kernelEndCol; kernelCol++)
            {
                sum += kernel[kernelRow, kernelCol] * img[imgRow - centerY + kernelRow, imgCol - centerX + kernelCol];
            }
        }
        response[imgRow, imgCol] = sum;
    }
}

A completely different approach

There remains a time complextiy of

$$O(n^2 \cdot k^2)$$

which quickly escalates execution time for large kernels, a situation we can't migitage by improvements on code level. To get around this bound, one can calculate the convolution in frequency domain. Consider this MATLAB code:

>> signal = [1, 2, 3, 4];
>> kernel = [.5, 1, .5];
>> conv(signal, kernel)
ans =
   0.50000   2.00000   4.00000   6.00000   5.50000   2.00000
>> ifft( fft([signal,0,0]) .* fft([kernel,0,0,0]) )
ans =
   0.50000   2.00000   4.00000   6.00000   5.50000   2.00000
>> ifft( fft([signal,0,0,0,0]) .* fft([kernel,0,0,0,0,0]) )
ans =
   0.50000   2.00000   4.00000   6.00000   5.50000   2.00000   0.00000  0.00000

A convolution can always be calculated by zero-padding both signal and kernel to a length of at least

$$n+k-1$$

Then apply the Fast Fourier Transform to each of them, multiply the results elementwise and calculate the inverse FFT. This leaves us with the time complexity of the 2D-FFT:

$$O((n+k)^2~log~(n+k))$$

That's a major improvement over the naive algorithm, except for small kernel sizes, such as 5 x 5. If you plan to use larger kernels, I strongly recommend to switch algorithms.

But wait, there's more: You don't even need to calculate the Gabor filter in spacial domain and plug it into a FFT. It can be transformed analytically. In fact, it's easier to represent it in frequency domain: It's just a shifted gaussian function. Distance and direction from origin determine modulation frequency and orientation.

How to implement the filter bank?

Zero-pad the input image to at least imagesize + kernelsize - 1. Round to the next higher power of two.
Apply the two-dimensional Fast Fourier Transform to the image.
Generate a set of Gabor filters directly in frequncy domain. You may also generate just one filter and apply it to different regions of the source image. Remember: Shifting determines modulation frequency and orientation.
Multiply the filters with the image, resulting in a set of filterd images.
Apply the inverse FFT to each of the filterd images.
Extract the valid region to get images of the original size. In the MATLAB example it's [2, 4, 6, 5.5].

"In addition, 14% of CPU ... ... ... by moving calculations and branches to higher level loops." --- You gave some sample code of loop-optimization. Where should I paste that code? I could't understand the location where I should place that code. — anonymous, Sep 1 at 7:32
"1. Zero-pad the input image to at least imagesize + kernelsize - 1. Round to the next higher power of two." --- If I only pad the input image, how can I multiply the kernel with that image? Wouldn't the image and the kernel be of different dimensions? — anonymous, Sep 1 at 7:37
"Generate a set of Gabor filters directly in frequency domain." --- Looks like a very good idea. Would you please let me know how can I generate that? — anonymous, Sep 1 at 7:45
@anonymous There's only one occation where you nest four loops, and thats where the code should go. It's in GaborFilter.Apply. — Rainer P., Sep 1 at 8:13
@anonymous The kernel must also be padded, but this is done implicitly. First, allocate a result array initalized to all zeros. Then, multiply some region of the source image with the kernel and write it to the corresponding region in the result image. Leave the rest of the result as is. This is equivalent to first pad the kernel and then multiply the full image. — Rainer P., Sep 1 at 8:17

asked	29 days ago
viewed	50 times
active	29 days ago

current community

your communities

more stack exchange communities

Gabor filter implementation

1 Answer 1

Speeding up your code

A completely different approach

How to implement the filter bank?

Your Answer

Not the answer you're looking for? Browse other questions tagged c# performance or ask your own question.

Hot Network Questions

current community

your communities

more stack exchange communities

Gabor filter implementation

1 Answer 1

Speeding up your code

A completely different approach

How to implement the filter bank?

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged c# performance or ask your own question.

Related

Hot Network Questions