Gabor滤波器

滤波器学习2014-01-0321:1910135人阅读评论(12)收藏举报图像处理gabor滤波本文的目的是用C实现生成Gabor模版，并对图像卷积。并简单提一下，Gabor滤波器在纹理特征提取上的应用。一、什么是Gabor函数（以下内容含部分翻译自维基百科）在图像处理中，Gabor函数是一个用于边缘提取的线性滤波器。Gabor滤波器的频率和方向表达同人类视觉系统类似。研究发现，Gabor滤波器十分适合纹理表达和分离。在空间域中，一个二维Gabor滤波器是一个由正弦平面波调制的高斯核函数。还有，生物学实验发现，Gabor滤波器可以很好地近似单细胞的感受野函数（光强刺激下的传递函数），什么视皮层内的超柱，bla...bla，总之是这方面仿生的数学模型。另外，网上有一种说法，gabor分为实部和虚部，用实部进行滤波后图像会平滑；虚部滤波后用来检测边缘。【来自百度知道某个大神的回答】，我查了文献，发现的确有人用Gabor的奇函数部分做边缘提取（《基于Gabor滤波器的边缘检测算法》无线电工程2000年第3卷第30期）。另外，从我的实验结果也有类似的发现。暂且认为这个对的吧。Gabor滤波器的脉冲响应，可以定义为一个正弦波（对于二维Gabor滤波器是正弦平面波）乘以高斯函数。由于乘法卷积性质，Gabor滤波器的脉冲响应的傅立叶变换是其调和函数的傅立叶变换和高斯函数傅立叶变换的卷积。该滤波器由实部和虚部组成，二者相互正交。一组不同频率不同方向的Gabor函数数组对于图像特征提取非常有用。下面给出二维Gabor函数的数学表达：复数表达：实数部分：虚数部分：其中：和下面介绍公式中各个参数的含义，及参数如何配置问题【都从老外那翻译来的】：波长（λ）：它的值以像素为单位指定，通常大于等于2.但不能大于输入图像尺寸的五分之一。方向（θ）：这个参数指定了Gabor函数并行条纹的方向，它的取值为0到360度相位偏移（φ）：它的取值范围为-180度到180度。其中，0he180度分别对应中心对称的center-on函数和center-off函数，而-90度和90度对应反对称函数。长宽比（γ）：空间纵横比，决定了Gabor函数形状（support，我翻译为形状）的椭圆率（ellipticity）。当γ=1时，形状是圆的。当γ1时，形状随着平行条纹方向而拉长。通常该值为0.5带宽（b）：Gabor滤波器的半响应空间频率带宽b和σ/λ的比率有关，其中σ表示Gabor函数的高斯因子的标准差，如下：σ的值不能直接设置，它仅随着带宽b变化。带宽值必须是正实数，通常为1，此时，标准差和波长的关系为：σ=0.56λ。带宽越小，标准差越大，Gabor形状越大，可见平行兴奋和抑制区条纹数量越多。下面给出，不同参数配置下的Gabor核函数效果图，大小均100*100：a.波长对比组【方向为：0，相位偏移量为：0，纵横比率为：0.5，带宽为：1，下图波长分别为5，10，15】b.方向对比组【波长为：10，相位偏移量为：0，空间纵横比为：0.5，带宽为：1，方向分别为：0，45，90】c.相位偏移量对比组【波长为：10，方向为：0，空间纵横比：0.5，带宽：1，相位偏移量分别为：0，180，-90，90】d.空间纵横比对比组【波长：10，相位偏移量：0，方向：0，带宽：1，空间纵横比分别为：0.5，1】e.带宽对比组【波长：10，方向：0，相位偏移量：0，空间纵横比：0.5，带宽分别为：0.5，1，2】二、gabor函数实现：matlab版本，我从pudn上找来的，但他的gabor函数，我没怎么看明白：gabor函数：[plain]viewplaincopy1.functiongabor_k=compute(x,y,f0,theta)2.r=1;g=1;3.x1=x*cos(theta)+y*sin(theta);4.y1=-x*sin(theta)+y*cos(theta);5.gabor_k=f0^2/(pi*r*g)*exp(-(f0^2*x1^2/r^2+f0^2*y1^2/g^2))*exp(i*2*pi*f0*x1);[plain]viewplaincopy1.%绘制一个Gabor滤波器的空域和频域函数图2.clear;3.x=0;4.theta=0;5.f0=0.2;6.fori=linspace(-15,15,50)7.x=x+1;8.y=0;9.forj=linspace(-15,15,50)10.y=y+1;11.z(y,x)=compute(i,j,f0,theta);12.end13.end14.x=linspace(-15,15,50);15.y=linspace(-15,15,50);16.surf(x,y,real(z))17.title('Gaborfilter:realcomponent');18.xlabel('x');19.ylabel('y');20.zlabel('z');21.figure(2);22.surf(x,y,imag(z))23.title('Gaborfilter:imaginarycomponent');24.xlabel('x');25.ylabel('y');26.zlabel('z');27.28.Z=fft2(z);29.u=linspace(-0.5,0.5,50);30.v=linspace(-0.5,0.5,50);31.figure(3);32.surf(u,v,abs(fftshift(Z)))33.title('Gaborfilter:frequencycomponent');34.xlabel('u');35.ylabel('v');36.zlabel('Z');运行结果：[plain]viewplaincopy1.%4个方向的Gabo滤波器通过图像显示2.clear;3.x=0;4.theta=pi*3/4;%用弧度0,pi/4,pi/2,pi*3/45.f0=0.2;6.fori=linspace(-15,15,50)7.x=x+1;8.y=0;9.forj=linspace(-15,15,50)10.y=y+1;11.z(y,x)=compute(i,j,f0,theta);12.end13.end14.z_real=real(z);15.m=min(z_real(:));16.z_real=z_real+abs(m);17.M=max(z_real(:));18.imshow(1/M*z_real);19.figure(2)20.z_imag=imag(z);21.m=min(z_imag(:));22.z_imag=z_imag+abs(m);23.M=max(z_imag(:));24.imshow(1/M*z_imag);运行效果：实数部分：虚数部分：[plain]viewplaincopy1.%4个方向的Gabor滤波器对lena进行滤波2.clear;3.I=imread('.\pic\lena.bmp');4.f0=0.2;5.count=0;6.fortheta=[0,pi/4,pi/2,pi*3/4];%用弧度0,pi/4,pi/2,pi*3/47.count=count+1;8.x=0;9.fori=linspace(-8,8,11)10.x=x+1;11.y=0;12.forj=linspace(-8,8,11)13.y=y+1;14.z(y,x)=compute(i,j,f0,theta);15.end16.end17.figure(count);18.filtered=filter2(z,I,'valid');19.f=abs(filtered);20.imshow(f/max(f(:)))21.end运行效果：好吧，不管他了。大概感受一下吧。由于我没看明白他的gabor函数怎么定义的，参数设置也不一样，实验结果很不相同，我希望我是对的，天地良心呐！！我只能按照维基百科给出的函数，编写了以下C代码：[cpp]viewplaincopy1.//my_gabor.cpp:定义控制台应用程序的入口点。2.//3.4.#includestdafx.h5.#includeiostream6.#includeopencv2/core/core.hpp7.#includeopencv2/highgui/highgui.hpp8.#includeopencv2/imgproc/imgproc.hpp9.#includemath.h10.#definePI3.141592611.#defineN412.usingnamespacestd;13.usingnamespacecv;14.15.voidm_filer(double*src,intheight,intwidth,double*mask_rel,double*mask_img,intmW,intmH,intk)16.{17.IplImage*tmp;18.doublea,b,c;19.charres[20];//保存的图像名称20.21.tmp=cvCreateImage(cvSize(width,height),IPL_DEPTH_8U,1);22.23.for(inti=0;iheight;i++)24.{25.for(intj=0;jwidth;j++)26.{27.a=0.0;28.b=0.0;29.c=0.0;30.//去掉靠近边界的行列，无法滤波，超出范围31.if(iint(mH/2)&&iheight-int(mH/2)&&jint(mW)&&jwidth-int(mW/2))32.{33.for(intm=0;mmH;m++)34.{35.for(intn=0;nmW;n++)36.{37.//printf(%f\n,src[(i+m-int(mH/2))*width+(j+n-int(mW))]);38.a+=src[(i+m-int(mH/2))*width+(j+n-int(mW))]*mask_rel[m*mW+n];39.b+=src[(i+m-int(mH/2))*width+(j+n-int(mW))]*mask_img[m*mW+n];40.//printf(%f,%f\n,a,b);41.}42.}43.}44.c=sqrt(a*a+b*b);45.c/=mW*mH;46.tmp-imageData[i*width+j]=(unsignedchar)c;47.}48.}49.sprintf(res,result%d.jpg,k);50.cvSaveImage(res,tmp);51.cvReleaseImage(&tmp);52.}53.54.int_tmain(intargc,_TCHAR*argv[])55.{56.IplImage*src;57.double*rel,*img,*src_data,xtmp,ytmp,tmp1,tmp2,tmp3,re,im;58.doubleTheta,sigma,Gamma,Lambda,Phi;//公式中的5个参数59.intgabor_height,gabor_width,x,y;60.61.src=cvLoadImage(test.jpg,CV_LOAD_IMAGE_GRAYSCALE);62.gabor_height=10;63.gabor_width=10;64.Gamma=1.0;65.Lambda=10.0;66.sigma=100;67.Phi=0;68.69.rel=(double*)malloc(sizeof(double)*gabor_width*gabor_height);//实数部分70.img=(double*)malloc(sizeof(dou