forstner算子提取特征点(matlab版)

clear;
clc;
[FileName,PathName,FilterIndex] = uigetfile({'*.tif;*.bmp;*.jpg;*.png','All image files'},'Select the input image');%获取所选图像的信息（文件名、路径等）
if FilterIndex
pic_orig=imread([PathName,FileName]); %读取原始标准图像
else
return
end
figure(1);
imshow(pic_orig);%显示原始标准图像
title('原始标准图像');

[x,y,z] = size(pic_orig); % 判断该图像是彩色图像还是灰度图像，如果是彩色图像把它转换为灰度图像
if(z~=1)
pic_gray = rgb2gray(pic_orig);
end
figure(2);
imshow(pic_gray); %显示灰度图像
title('灰度图像');

% level=graythresh(pic_gray);%二值图像的阈值
pic_bw=im2bw(pic_gray,0.9); %把灰度图像转换成二值图像,阈值可调
figure(3);
imshow(pic_bw);
title('二值化图像');

%第一步：先用4个方向的差分算子提取初选点
result=zeros(x,y);%特征提取结果
for i=2:x-1
for j=2:y-1
dg1=abs(pic_bw(i,j)-pic_bw(i+1,j));
dg2=abs(pic_bw(i,j)-pic_bw(i,j+1));
dg3=abs(pic_bw(i,j)-pic_bw(i-1,j));
dg4=abs(pic_bw(i,j)-pic_bw(i,j-1));
dg=[dg1,dg2,dg3,dg4];
temp=sort(dg); %对4个差分算子进行升序排列
if temp(3)==1 %有任意两个大于阈值，则该像素有可能是一特征点
result(i,j)=255; %是初选点
else
result(i,j)=0; %不是初选点
end
end
end
figure(4);
imshow(result);
title('初选点');

%第二步：在以初选点为中心的3*3的窗口中计算协方差矩阵与圆度
wMatrix=zeros(x,y); %权重矩阵
Tq=0.8; %阈值，可设置
for i=2:x-1
for j=2:y-1
if result(i,j)==255 %如果是初选点
gu2=0.0;
gv2=0.0;
guv=0.0;
for ii=i-1:i
for jj=j-1:j
gu2=gu2+(pic_bw(ii+1,jj+1)-pic_bw(ii,jj))^2;
gv2=gv2+(pic_bw(ii,jj+1)-pic_bw(ii+1,jj))^2;
guv=guv+(pic_bw(ii+1,jj+1)-pic_bw(ii,jj))*(pic_bw(ii,jj+1)-pic_bw(ii+1,jj));
end
end
DetN=gu2*gv2-guv^2;
trN=gu2+gv2;
q=4*DetN/(trN*trN);
%第三步：设定阈值Tq,若满足则计算权值
if q > Tq
wMatrix(i,j)=DetN/trN;
else
result(i,j)=0;
end
end
end
end
figure(5);
imshow(result);
title('进一步提取的特征点');

%第四步：以权值为基础，在一定窗口内抑制局部非最大值候选点,取出局部极大值点
vwsize=5; %选择5*5的窗口
wradius=floor(vwsize/2);
for i=wradius+1:x-wradius
for j=wradius+1:y-wradius
if result(i,j)==255
tempiv=wMatrix(i-wradius:i+wradius,j-wradius:j+wradius);
tempsort=sort(tempiv(:),'descend');%将区域内像素按从大至小排列
if wMatrix(i,j)==tempsort(1) && wMatrix(i,j)~=tempsort(2)%排除整个区域权值相等的情况
;

else
result(i,j)=0;
end
end
end
end
figure(6);
imshow(result);
title('forstner提取的最终匹配点');