5
$\begingroup$

For the application of Gabor filters in feature-based (e.g. lines and edges in images) digital image processing, it is desired that the DC component is 0. Why is this so? What effect does it have mathematically if the DC component of the filter is not 0? Thanks in advance!

$\endgroup$
5
$\begingroup$

In short: If a filter have DC component it means that for a constant region of the image, the filter response will depend on the grey value of that region. Constant region with higher values ==> higher response.

One usually wants higher responses only in case of edges/corners and not in the case of bright but constant regions.

As a picture worth 1000 words, here is an example. I took a simple image consisting of two horisontal lines, each one in different background. I took the Gabor and LogGabor filter. LogGabor is a variant that have 0 DC response.

I filtered the image with a Gabor filter (tuned to respond to horizontal signals) and then, the same image, with LogGabor filter (tuned with exactly the same parameters as Gabor filter).

Original image:

Original image

Gabor convoluted: gabor

LogGabor convoluted: logGabor

Another example of Gabor filter, with a bandwidth of 0.7 octaves. gabor0.7

Now, let's ignore the edges for a moment and focus on the constant regions. In the second image, you see different responses for the constant regions! Darker region have a lower response than the lighter one. For the LogGabor case (third image) you see that the response is the same regardless of the original intensity.

Note: The response values were scaled in 0-255 interval, so, neutral grey means 0 response, darker means negative values and lighter positive ones. The bandwidth is set to 2 octaves (except for the 4'th image)

EDIT

According to a paper by Boukerroui et al., for bandwidths around 0.7 the DC response is close to the computational errors ($10^{-6}$). The 4'th figure demonstrate this.

My personal advice to you is to use another bandpass filter. I strongly recommend Gauss Derivatives or Difference of Gaussians (in this order) if it is necessary for you to stay in spatial domain. If you don't have this restriction, I personally prefer LogGabor.

From my experience there is no net advantage or disadvantage in terms of end result between logGabor, Gaussian Derivatives, Difference of Gaussians or more exotic Cauchy filters. BUT this might vary, depending on your practical application.

As for the code, here is a snippet of my Gabor/LogGabor implementation. The code is in Matlab, the convolutions are done in Fourier space. The complexity is $O(n log n)$ and it depends only on the image dimensions, not on the filter dimensions. The code follows the guidelines in Boukerroui et al.. It is developed by me and released under MIT licence (not shown here) so feel free to use it as you see fit.

Img must be a MxNx1 of doubles, (preferably between 0 and 1 but it should work with 0-255 too); P0 is the tuning period (wavelength), that is, $P_0=1/f$, where f is tuning frequency; Orient is the filter orientation, expressed in radians; FBW frequency bandwidth in octaves; ABW the angle bandwidth.

function [ ReConv,ImConv ] = GaborConv( Img, P0,Orient, FBW, ABW,varargin )
 ReConv=Img; ImConv=Img;
 ImgFFt=fft2(Img);
 lg2pi=sqrt(log(2)/pi);
 A=P0*lg2pi*(2^FBW+1)/(2^FBW-1); B=P0*lg2pi*(1/tan(0.5*ABW));
 K=1/(A*B); KF=2*pi/P0;
 [N,M]=size(Img); u0=1/P0*cos(Orient); v0=1/P0*sin(Orient);
 FFilt=zeros([M,N]);
 stx=1/(M-1); sty=1/(N-1);
 [U V]=meshgrid(-0.5:stx:0.5,-0.5:sty:0.5);
 Ur=(U-u0)*cos(Orient)+(V-v0)*sin(Orient);
 Vr=-(U-u0)*sin(Orient)+(V-v0)*cos(Orient);
 FFilt=exp(-pi*(Ur.*Ur.*A.*A+Vr.*Vr.*B.*B)) ;
 FFilt=ifftshift(FFilt);
 normFact=sum(sum((abs(ifft2(FFilt)))));
 RezFreq=FFilt.*ImgFFt/normFact; RezSpace=ifft2(RezFreq);
 ReConv=real(RezSpace); ImConv=imag(RezSpace);
end

Log Gabor method:

function [  ReConv,ImConv ] = LogGaborConv( Img, P0,Orient, FBW, ABW,useCos ).
 ReConv=Img;
 ImConv=Img;
 ImgFFt=fft2(Img);
 %Compute the k_beta according to the desired frequency bandwidth
 kBeta=exp(-1/4*sqrt(2*log(2))*FBW); 
 gaussSD=sqrt(ABW^2/(2*log(2))); %compute Gauss SD to ensure that Gauss(ABW,0,SD)==0.5
 lnkB2=2*log(kBeta)^2;
 nc=exp(-1/8*lnkB2/2)*sqrt(-2*sqrt(pi)/(1/P0*log(kBeta)));
 [N M]=size(Img);
 stx=1/(M-1);
 sty=1/(N-1);
 [Ug Vg]=meshgrid(-0.5:stx:0.5,-0.5:sty:0.5);
 Radius=sqrt(Ug.^2+Vg.^2);
 Radius = ifftshift(Radius);
 Radius(1,1)=1;
 Radius = fftshift(Radius);
 %prepare a low pass filter, as proposed by Kovesi
 lowpass=1./(1+(Radius./0.45).^(30));
 %Compute the bandpass filter
 logRad=nc*exp(-log(Radius*P0).^2/lnkB2);
 RadialBP=logRad.*lowpass;
 %Compute the angle component
 Theta=atan2(Vg,Ug);
 sintheta = sin(Theta);
 costheta = cos(Theta);
 ds = sintheta * cos(Orient) - costheta * sin(Orient);
 dc = costheta * cos(Orient) + sintheta * sin(Orient);
 dtheta = abs(atan2(ds,dc));
 if ~exist('useCos','var')
     useCos=0;
 end
 if useCos==0
     %Use a gauss angle
     AngFilter=exp(-dtheta.^2/(2*gaussSD^2));
 else
     %Use a cosine angle
     AngFilter=cos(dtheta).^2.*heaviside(dtheta./pi+0.5).*heaviside(0.5-dtheta./pi);
 end
 logGaborFilt=RadialBP.*AngFilter;
 logGaborFilt=ifftshift(logGaborFilt);
 logGaborFilt(1,1)=0;
 normFact=sum(sum((abs(ifft2(logGaborFilt)))));
 RezFreq=logGaborFilt.*ImgFFt/normFact;
 RezSpace=ifft2(RezFreq);
 ReConv=real(RezSpace);
 ImConv=imag(RezSpace);
end

As for how the DC response varies with various parameters and why, a comprehensive discussion is in Boukerroui et al., equations (19) and (20) at pages 75, 76. Take a peek at Appendix C (pag 101). I would like to redirect you to the paper rather than just copy/paste the relevant paragraphs.

Hope it helps!

$\endgroup$
  • $\begingroup$ Hi @visoft, thanks for your answer! anyway is it possible to make the gabor filter 0 DC to improve its performance? $\endgroup$ – freak_warrior Feb 4 '14 at 13:21
  • 1
    $\begingroup$ In theory, if you make the frequency bandwidth small enough (~0.7) the DC response will be low (close to the computational errors). I don't use it because lower bandwidth means larger spatial support (larger filter) I'll edit the post later with an example and some references. $\endgroup$ – visoft Feb 4 '14 at 13:28
  • $\begingroup$ Can you show mathematically why does the DC response varies with the values of the image if it is nonzero? Also can you show me your code for Gabor and log Gabor respectively?:) $\endgroup$ – freak_warrior Feb 4 '14 at 15:07
  • $\begingroup$ @freak_warrior I added more details in the answer. $\endgroup$ – visoft Feb 4 '14 at 16:31
  • $\begingroup$ p.s. A good piece of code is in your previous question: dsp.stackexchange.com/questions/13907/… $\endgroup$ – visoft Feb 4 '14 at 16:34
4
$\begingroup$

I would like to add some notes to @visoft's answer.

Log-gabor filter is a very good one alternative to gabor filter has no DC component. Yet in ordinary gabor filter it is also possible to remove the DC component. See the Matlab code below:

function [GABOUT1,GABOUT2]=gaborfilter(I,S,F,W,P)

if isa(I,'double')~=1
    I=double(I);
end
size=fix(1.5/S); % exp(-1.5^2*pi) < 0.1%

k=1;
for x=-size:size
    for y=-size:size
        G1(size+x+1,size+y+1)=k*exp(-pi*S^2*(x*x+y*y))*...
            (exp(sqrt(-1)*(2*pi*F*(x*cos(W)+y*sin(W))+P)));
        G2(size+x+1,size+y+1)=k*exp(-pi*S^2*(x*x+y*y))*...
            (exp(sqrt(-1)*(2*pi*F*(x*cos(W)+y*sin(W))+P))-exp(-pi*(F/S)^2+sqrt(-1)*P));
    end
end

GABOUT1=conv2(I,double(G1),'same');
GABOUT2=conv2(I,double(G2),'same');

In the code exp(-pi*(F/S)^2+sqrt(-1)*P)) is the DC component due to the overlapping. And GABOUT1 is the filtered image with DC component, while GABOUT2 is the filtered image with the compensation of DC component.

The issue with DC component is that one cannot achieve both broad spectral information and maximal spatial localization. As a result, you may over-estimate the low-frequency components but under-estimate the high-frequency components in a image. And this is especially reflected if you observe the magnitude of the filtered image. But there is no evident difference if you are observing the imag part of the image.

Check one example below:

I=imread('cat.jpg');
I=rgb2gray(I);
imshow(I),figure,
[G1,G2]=gaborfilter(I,0.5,0.1,0,0);
I1=abs(G1);
I1=I1/max(I1(:))*256;
imshow(uint8(I1))

I2=abs(G2);
I2=I2/max(I2(:))*256;
figure,imshow(uint8(I2))

Compared with the original image (left), the filtered one with DC component (middle) deteriorates at the high frequency, while the one without DC component (right) enhances the high-frequency details (tiny edges at various orientations) very much.

You can observe the imag part of the filtered image by changing the code to I1=imag(G1); and I2=imag(G2); and you may not see that much difference, i.e., the influence of DC component is not that much at the imag part of the filtered image.

$\endgroup$
  • $\begingroup$ Hi lennon310: thanks for your answer. BTW u posted an answer for the prediction of fine-level wavelet coefficients. Could you explain more for the derivation of the formula for that question? $\endgroup$ – freak_warrior Feb 5 '14 at 1:23
  • $\begingroup$ Thank you freak_warrior. I couldn't find that post, seems like being deleted for some reason. Actually that formula is from a publication, and I noticed later that they made some assumptions on that, so probably that's not the coefficient for general case. Sorry about that. $\endgroup$ – lennon310 Feb 5 '14 at 1:30
  • 1
    $\begingroup$ @lennon310 The DC removal solution looks like Ronse et al. proposal. Am I right? $\endgroup$ – visoft Feb 7 '14 at 15:41
  • 1
    $\begingroup$ @freak_warrior I see that nobody answered your remark for past 2 days, so here is a quick and dirty attempt: In the Appendix C of the reference paper I put in my answer, there are two methods that deal with DC response. Both are fairly complex and don't not completely solves the problem. IMHO whatever you do in the Fourier space or to the spatial filter itself will alter its properties in an unpredictable way and you will no longer have a Gabor filter. $\endgroup$ – visoft Feb 7 '14 at 15:43
  • 1
    $\begingroup$ @visoft Movellan, J. R. - Tutorial on Gabor Filters. Tech. rep., 2002. $\endgroup$ – lennon310 Feb 7 '14 at 15:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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