How to apply Gabor wavelets to an image?

A typical use of Gabor filters is to calculate filter responses at each of several orientations, e.g. for edge detection.

You can convolve a filter with an image using the Convolution Theorem, by taking the inverse Fourier transform of the element-wise product of the Fourier transforms of the image and the filter. Here is the basic formula:

%# Our image needs to be 2D (grayscale)
if ndims(img) > 2;
    img = rgb2gray(img);
%# It is also best if the image has double precision
img = im2double(img);

[m,n] = size(img);
[mf,nf] = size(GW);
GW = padarray(GW,[n-nf m-mf]/2);
GW = ifftshift(GW);
imgf = ifft2( fft2(img) .* GW );

Typically, the FFT convolution is superior for kernels of size > 20. For details, I recommend Numerical Recipes in C, which has a good, language agnostic description of the method and its caveats.

Your kernels are already large, but with the FFT method they can be as large as the image, since they are padded to that size regardless. Due to the periodicity of the FFT, the method performs circular convolution. That means that the filter will wrap around the image borders, so we have to pad the image itself as well to eliminate this edge effect. Finally, since we want the total response to all the filters (at least in a typical implementation), we need to apply each to the image in turn, and sum the responses. Usually one uses just 3 to 6 orientations, but it is also common to do the filtering at several scales (different kernel sizes) as well, so in that context a larger number of filters is used.

You can do the whole thing with code like this:

img = im2double(rgb2gray(img)); %# 
[m,n] = size(img); %# Store the original size.

%# It is best if the filter size is odd, so it has a discrete center.
R = 127; C = 127;

%# The minimum amount of padding is just "one side" of the filter.
%# We add 1 if the image size is odd.
%# This assumes the filter size is odd.
pR = (R-1)/2;
pC = (C-1)/2;
if rem(m,2) ~= 0; pR = pR + 1; end;
if rem(n,2) ~= 0; pC = pC + 1; end;
img = padarray(img,[pR pC],'pre'); %# Pad image to handle circular convolution.

GW = {}; %# First, construct the filter bank.
for v = 0 : 4
    for u = 1 : 8
        GW =  [GW {GaborWavelet(R, C, Kmax, f, u, v, Delt2)}];

%# Pad all the filters to size of padded image.
%# We made sure padsize will only be even, so we can divide by 2.
padsize = size(img) - [R C];
GW = cellfun( ...
        @(x) padarray(x,padsize/2), ...
        GW, ...

imgFFT = fft2(img); %# Pre-calculate image FFT.

for i=1:length(GW)
    filter = fft2( ifftshift( GW{i} ) ); %# See Numerical Recipes.
    imgfilt{i} = ifft2( imgFFT .* filter ); %# Apply Convolution Theorem.

%# Sum the responses to each filter. Do it in the above loop to save some space.
imgS = zeros(m,n);
for i=1:length(imgfilt)
    imgS = imgS + imgfilt{i}(pR+1:end,pC+1:end); %# Just use the valid part.

%# Look at the result.

Keep in mind that this is essentially the bare minimum implementation. You can choose to pad with replicates of the border instead of zeros, apply a windowing function to the image or make the pad size larger in order to gain frequency resolution. Each of these is a standard augmentation to the technique I’ve outlined above and should be trivial to research via Google and Wikipedia. Also note that I haven’t added any basic MATLAB optimizations such as pre-allocation, etc.

As a final note, you may want to skip the image padding (i.e. use first code example) if your filters are always much smaller than the image. This is because, adding zeros to an image creates an artificial edge feature where the padding begins. If the filter is small, the wraparound from circular convolution doesn’t cause problems, because only the zeros in the filter’s padding will be involved. But as soon as the filter is large enough, the wraparound effect will become severe. If you have to use large filters, you may have to use a more complicated padding scheme, or crop the edge of the image.

Leave a Comment