1 简介
D.Gabor 1946年提出
窗口Fourier变换，为了由信号的Fourier变换提取局部信息，引入了时间局部化的窗函数。
由于窗口Fourier变换只依赖于部分时间的信号，所以，现在窗口Fourier变换又称为短时Fourier变换，这个变换又称为Gabor变换。
1) Gabor优点
Gabor小波与人类视觉系统中简单细胞的视觉刺激响应非常相似。它在提取目标的局部空间和频率域信息方面具有良好的特性。虽然Gabor小波本身并不能构成正交基，但在特定参数下可构成紧框架。Gabor小波对于图像的边缘敏感，能够提供良好的方向选择和尺度选择特性，而且对于光照变化不敏感,能够提供对光照变化良好的适应性。上述特点使Gabor小波被广泛应用于视觉信息理解。
Gabor滤波器和脊椎动物视觉皮层感受野响应的比较：第一行代表脊椎动物的视觉皮层感受野，第二行是Gabor滤波器，第三行是两者的残差。可见两者相差极小。Gabor滤波器的这一性质，使得其在视觉领域中经常被用来作图像的预处理。
2) Gabor定义
① 具体窗函数――Gaussaion的 Gabor变换定义式
Gabor变换的基本思想：把信号划分成许多小的时间间隔，用傅里叶变换分析每一个时间间隔，以便确定信号在该时间间隔存在的频率。其处理方法是对f(t)加一个滑动窗，再作傅里叶变换。
② 窗口的宽高关系
经理论推导可以得出：高斯窗函数条件下的窗口宽度与高度，且积为一固定值。
3) 离散Gabor变换的一般求法
① 首先选取核函数
可根据实际需要选取适当的核函数。如，如高斯窗函数；
② 离散Gabor变换的表达式
4) Gabor变换的解析理论
对偶函数可以使计算更为简洁方便。
5) 适用条件
① 临界采样Gabor展开要求条件：TΩ=2π；
② 过采样展开要求条件：TΩ≤2π；
当TΩ＞2π时，欠采样Gabor展开，已证明会导致数值上的不稳定。
6) 应用
① 暂态信号检测
如果对信号波形有一定的先验知识且可以据此选取合适的基函数，可以用Gabor变换对信号作精确的检测统计计量。
② 图象分析与压缩
二维Gabor变换可以应用到图象分析与压缩中。
3. 二维Gabor滤波器
用Gabor 函数形成的二维Gabor 滤波器具有在空间域和频率域同时取得最优局部化的特性，因此能够很好地描述对应于空间频率(尺度)、空间位置及方向选择性的局部结构信息。
Gabor滤波器的频率和方向表示接近人类视觉系统对于频率和方向的表示，并且它们常备用于纹理表示和描述。
在图像处理领域，Gabor滤波器是一个用于边缘检测的线性滤波器。
在空域，一个2维的Gabor滤波器是一个正弦平面波和高斯核函数的乘积。
Gabor滤波器是自相似的，也就是说，所有Gabor滤波器都可以从一个母小波经过膨胀和旋转产生。
实际应用中，Gabor滤波器可以在频域的不同尺度，不同方向上提取相关特征。
Gabor滤波器的傅里叶变换：峰值响应在复正弦的空域频率(u0,v0)：
Gabor滤波器示意图，3种角度5种方向：
2 部分代码
b1=[0 1 0
1 1 1
0 1 0];
b2=[0 0 1 0 0
0 1 1 1 0
1 1 1 1 1
0 1 1 1 0
0 0 1 0 0];
b31=[0 0 0
1 1 1
0 0 0];
b32=[0 0 1
0 1 0
1 0 0];
b33=[0 1 0
0 1 0
0 1 0];
b34=[1 0 0
0 1 0
0 0 1];
I = imread('C:\Documents and Settings\Administrator\桌面\matlab code\邵海峰\matlab代码\原始指纹图像2.jpg');
I=rgb2gray(I);
I=double(I);
figure,imshow(I,[])
xs1=size(I,1);
xs2=size(I,2);
K1=sum(sum(I));
M1=K1/(xs1*xs2);
M2=ones(xs1,xs2);
M0=M2*M1;
M3=zeros(xs1,xs2);
for i=1:xs1
for j=1:xs2
M3(i,j)=(I(i,j)-M0(i,j))^2; %求一个平均值
end
end
K2=sum(sum(M3));
V0=K2/(xs1*xs2);
P1=ones(fix(xs1/15),fix(xs2/15));
P2=ones(15,15);
P3=ones(fix(xs1/15),fix(xs2/15));
for i=1:xs1
for j=1:xs2
a=mod(i,15);
b=mod(j,15);
if a==0
a=15;
end
if b==0
b=15;
end
P2(a,b)=I(i,j);
if a==15&&b==15
PK1=sum(sum(P2));
PM1=PK1/(15*15);
PM2=ones(15,15);
PM0=PM2*PM1;
PM3=zeros(15,15);
for m=1:15
for n=1:15
PM3(m,n)=(PM2(m,n)-PM0(m,n))^2;
end
end
PK2=sum(sum(PM3));
PV0=PK2/(15*15);
P1(fix(i/15),fix(j/15))=PM1;
P3(fix(i/15),fix(j/15))=PV0;
end
end
end
I_n=zeros(xs1,xs2);
for i=1:xs1
for j=1:xs2
m=ceil(i/15);
n=ceil(j/15);
if I(i,j)>M0(i,j)&&m~=ceil(xs1/15)&&n~=ceil(xs2/15)
I_n(i,j)=M0(i,j)+sqrt(V0*(I(i,j)-P1(m,n))^2/P3(m,n));
end
if I(i,j)<=M0(i,j)&&m~=ceil(xs1/15)&&n~=ceil(xs2/15)
I_n(i,j)=M0(i,j)-sqrt(V0*(P1(m,n)-I(i,j))^2/P3(m,n));
end
end
end
bw1= I;
vertical=[ 1 2 1
0 0 0
-1 -2 -1];
horizontal=[ -1 0 1
-2 0 2
-1 0 1];
fx=conv2(bw1,vertical,'same');
fy=conv2(bw1,horizontal,'same');
mubanver=15;mubanhor=15; %模板的垂直和水平大小
middle_mu1=size(bw1,1)/mubanver;
middle_mu2=size(bw1,2)/mubanhor;
middle_mu1=fix(middle_mu1);
middle_mu2=fix(middle_mu2);
thea= zeros(middle_mu1, middle_mu2);
thea1= zeros(middle_mu1, middle_mu2);
sum1= zeros(middle_mu1, middle_mu2);
sum2= zeros(middle_mu1, middle_mu2);
sum3= zeros(middle_mu1, middle_mu2);
for i=1:middle_mu1
for j=1:middle_mu2
for y=1:mubanver
for x=1:mubanhor
sum1(i,j)=sum1(i,j)+2*fx(y+(i-1)*(mubanver),x+(j-1)*(mubanhor))*fy(y+(i-1)*(mubanver),x+(j-1)*(mubanhor));
%2Gxy
sum2(i,j)=sum2(i,j)+fx(y+(i-1)*(mubanver),x+(j-1)*(mubanhor))*fx(y+(i-1)*(mubanver),x+(j-1)*(mubanhor)) - fy(y+(i-1)*(mubanver),x+(j-1)*(mubanhor))*fy(y+(i-1)*(mubanver),x+(j-1)*(mubanhor));
%Gxx-Gyy
sum3(i,j)=sum3(i,j)+fx(y+(i-1)*(mubanver),x+(j-1)*(mubanhor))*fx(y+(i-1)*(mubanver),x+(j-1)*(mubanhor)) + fy(y+(i-1)*(mubanver),x+(j-1)*(mubanhor))*fy(y+(i-1)*(mubanver),x+(j-1)*(mubanhor));
%Gxx+Gyy
end
end
if sum2(i,j)>=0
fi=1/2*atan(sum1(i,j)/sum2(i,j));
end
if sum2(i,j)<0 && sum1(i,j)>=0
fi=1/2*(atan(sum1(i,j)/sum2(i,j))+pi);
end
if sum2(i,j)<0 && sum1(i,j)<0
fi=1/2*(atan(sum1(i,j)/sum2(i,j))-pi);
end
thea(i,j)=pi-fi; %方向角度
end
end
for i=1:size(sum3,1)
for j=1:size(sum3,2)
if sum3(i,j)==0
coh(i,j)=1;
else
coh(i,j)=sqrt(sum2(i,j)*sum2(i,j)+sum1(i,j)*sum1(i,j))/sum3(i,j);
end
end
end
figure,imshow(coh,[]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%根据置信度改变方向%%%%%%%%%%%%%%%%%%%%%%%%
[X,Y] = meshgrid(size(thea,1):-1:1,1:1:size(thea,2));
xx=X';yy=Y';
for i=1:size(thea,1)
for j=1:size(thea,2)
sum2(i,j)=cos(thea(i,j));
sum1(i,j)=sin(thea(i,j));
end
end
figure,quiver(yy,xx,sum2,sum1);
xs1=size(I,1);
xs2=size(I,2);
thea_jz = thea;
I=double(I);
bw1= I;
mubanver=15;mubanhor=15; %模板的垂直和水平大小
middle_mu1=size(bw1,1)/mubanver;
middle_mu2=size(bw1,2)/mubanhor;
middle_mu1=fix(middle_mu1);
middle_mu2=fix(middle_mu2);
for i=1:middle_mu1
for j=1:middle_mu2
if thea_jz(i,j)>=7/16*pi&&thea_jz(i,j)<9/16*pi
thea_jz(i,j)=1/2*pi;
end
if thea_jz(i,j)>=9/16*pi&&thea_jz(i,j)<11/16*pi
thea_jz(i,j)=5/8*pi;
end
if thea_jz(i,j)>=11/16*pi&&thea_jz(i,j)<13/16*pi
thea_jz(i,j)=3/4*pi;
end
if thea_jz(i,j)>=13/16*pi&&thea_jz(i,j)<15/16*pi
thea_jz(i,j)=7/8*pi;
end
if thea_jz(i,j)>=15/16*pi&&thea_jz(i,j)<17/16*pi
thea_jz(i,j)=pi;
end
if thea_jz(i,j)>=17/16*pi&&thea_jz(i,j)<19/16*pi
thea_jz(i,j)=9/8*pi;
end
if thea_jz(i,j)>=19/16*pi&&thea_jz(i,j)<21/16*pi
thea_jz(i,j)=5/4*pi;
end
if thea_jz(i,j)>=21/16*pi&&thea_jz(i,j)<23/16*pi
thea_jz(i,j)=11/8*pi;
end
if thea_jz(i,j)>=23/16*pi&&thea_jz(i,j)<25/16*pi
thea_jz(i,j)=1/2*pi;
end
end
end
for i=1:middle_mu1
for j=1:middle_mu2
for a=1:15
for b=1:15
c=a+(i-1)*15;
d=b+(j-1)*15;
theta_all(c,d)=thea_jz(i,j);
end
end
end
end
%%%%%%%%%%%%%%%%%%%%%thea=1/2*pi
Sx1=5;
Sy1=5;
theta1=1/2*pi;
f1=0.125;
if isa(I,'double')~=1
I = double(I);
end
for x = -fix(Sx1):fix(Sx1)
for y = -fix(Sy1):fix(Sy1)
xPrime = x * cos(theta1) + y * sin(theta1);
yPrime = y * cos(theta1) - x * sin(theta1);
G1(fix(Sx1)+x+1,fix(Sy1)+y+1) = exp(-.5*((xPrime/Sx1)^2+(yPrime/Sy1)^2))*cos(2*pi*f1*xPrime);
end
end
Imgabout = conv2(I,double(imag(G1)),'same');
Regabout = conv2(I,double(real(G1)),'same');
Gabout1 = sqrt(Imgabout.*Imgabout + Regabout.*Regabout)/10;
for i=1:xs1-8
for j=1:xs2-1
if theta_all(i,j)~=1/2*pi&&theta_all(i,j)~=11/8*pi
Gabout1(i,j)=0;
end
if Gabout1(i,j)>160
Gabout1(i,j)=256;
end
end
end
% Gabout1=uint8(Gabout1);
% F1=im2bw(Gabout1,85/255);1
% h=ones(3,3)/15;
% Gabout1=imfilter(F1,h);
%figure,imshow(Gabout1,[]);
% %%%%%%%%%%%%%%%%%%%% thea=5/8*pi
Sx2=7;
Sy2=7;
theta2=5/8*pi;
f2=0.1;
if isa(I,'double')~=1
I = double(I);
end
for x = -fix(Sx2):fix(Sx2)
for y = -fix(Sy2):fix(Sy2)
xPrime = x * cos(theta2) + y * sin(theta2);
yPrime = y * cos(theta2) - x * sin(theta2);
G2(fix(Sx2)+x+1,fix(Sy2)+y+1) = exp(-.5*((xPrime/Sx2)^2+(yPrime/Sy2)^2))*cos(2*pi*f2*xPrime);
end
end
Imgabout = conv2(I,double(imag(G2)),'same');
Regabout = conv2(I,double(real(G2)),'same');
Gabout2 = sqrt(Imgabout.*Imgabout + Regabout.*Regabout)/10;
for i=1:xs1-8
for j=1:xs2-1
if theta_all(i,j)~=5/8*pi
Gabout2(i,j)=0;
end
if Gabout2(i,j)>220
Gabout2(i,j)=256;
end
end
end
% Gabout2=uint8(Gabout2);
% F1=im2bw(Gabout2,20/255);
% h=ones(3,3)/15;
% Gabout2=imfilter(F1,h);
%figure,imshow(Gabout2,[]);
%%%%%%%%%%%%%%%%%%%% thea=3/4*pi xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
% Sx3=13;
% Sy3=13;
% theta3=3/4*pi;
% f3=0.09;
% if isa(I,'double')~=1
% I = double(I);
% end
% for x = -fix(Sx3):fix(Sx3)
% for y = -fix(Sy3):fix(Sy3)
% xPrime = x * cos(theta3) + y * sin(theta3);
% yPrime = y * cos(theta3) - x * sin(theta3);
% G3(fix(Sx3)+x+1,fix(Sy3)+y+1) = exp(-.5*((xPrime/Sx3)^2+(yPrime/Sy3)^2))*cos(2*pi*f3*xPrime);
% end
% end
% Imgabout = conv2(I,double(imag(G3)),'same');
% Regabout = conv2(I,double(real(G3)),'same');
% Gabout3 = sqrt(Imgabout.*Imgabout + Regabout.*Regabout)/10;
% for i=1:xs1-1
% for j=1:xs2-1
% if theta_all(i,j)~=3/4*pi
% Gabout3(i,j)=0;
% end
% % if Gabout3(i,j)>150
% % Gabout3(i,j)=256;
% % end
% end
% end
% Gabout3=uint8(Gabout3);
% F1=im2bw(Gabout3,88/255);
% h=ones(3,3)/15;
% Gabout3=imfilter(F1,h);
% figure,imshow(Gabout3,[]);
%%%%%%%%%%%%%%%%%%%% thea=7/8*pi
Sx4=7;
Sy4=7;
theta4=7/8*pi;
f4=0.12;
if isa(I,'double')~=1
I = double(I);
end
for x = -fix(Sx4):fix(Sx4)
for y = -fix(Sy4):fix(Sy4)
xPrime = x * cos(theta4) + y * sin(theta4);
yPrime = y * cos(theta4) - x * sin(theta4);
G4(fix(Sx4)+x+1,fix(Sy4)+y+1) = exp(-.5*((xPrime/Sx4)^2+(yPrime/Sy4)^2))*cos(2*pi*f4*xPrime);
end
end
Imgabout = conv2(I,double(imag(G4)),'same');
Regabout = conv2(I,double(real(G4)),'same');
Gabout4 = sqrt(Imgabout.*Imgabout + Regabout.*Regabout)/10;
for i=1:xs1-8
for j=1:xs2-1
if theta_all(i,j)~=7/8*pi&&theta_all(i,j)~=3/4*pi&&theta_all(i,j)~=pi
Gabout4(i,j)=0;
end
if Gabout4(i,j)>150
Gabout4(i,j)=256;
end
end
end
% Gabout4=uint8(Gabout4);
% F1=im2bw(Gabout4,30/255);
% h=ones(3,3)/15;
% Gabout4=imfilter(F1,h);
%figure,imshow(Gabout4,[]);
%
% %%%%%%%%%%%%%%%%%%%% thea=pi xxxxxxxxxxxxxxxxxxxxxxxxxxx
% Sx5=6;
% Sy5=6;
% theta5=pi;
% f5=0.08;
% if isa(I,'double')~=1
% I = double(I);
% end
% for x = -fix(Sx5):fix(Sx5)
% for y = -fix(Sy5):fix(Sy5)
% xPrime = x * cos(theta5) + y * sin(theta5);
% yPrime = y * cos(theta5) - x * sin(theta5);
% G5(fix(Sx5)+x+1,fix(Sy5)+y+1) = exp(-.5*((xPrime/Sx5)^2+(yPrime/Sy5)^2))*cos(2*pi*f5*xPrime);
% end
% end
%
% Imgabout = conv2(I,double(imag(G5)),'same');
% Regabout = conv2(I,double(real(G5)),'same');
% Gabout5 = sqrt(Imgabout.*Imgabout + Regabout.*Regabout)/10;
%
% for i=1:xs1-1
% for j=1:xs2-1
% if theta_all(i,j)~=pi
% Gabout5(i,j)=0;
% end
% end
% end
% Gabout5=uint8(Gabout5);
% F1=im2bw(Gabout5,230/255);
% h=ones(3,3)/15;
% Gabout5=imfilter(F1,h);
% figure,imshow(Gabout5,[]);
for i=1:xs1
for j=1:xs2
gabout_18=[Gabout1(i,j);Gabout2(i,j);Gabout4(i,j);Gabout6(i,j)];
gabout_last(i,j)=max(gabout_18);
if gabout_last(i,j)~=256
gabout_last(i,j)=0;
end
if i<20||j<30||i>xs1-20||j>xs2-20
gabout_last(i,j)=256;
end
end
end
figure,imshow(gabout_last,[]);
3 仿真结果
4 参考文献
[1]贡玉南, 华建兴, 黄秀宝. 基于匹配GabOr滤波器的规则纹理缺陷检测方法[J]. 中国图象图形学报, 2001, 006(007):624-628.
博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真，相关matlab代码问题可私信交流。
部分理论引用网络文献，若有侵权联系博主删除。