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

clear;clc;[FileName,PathName,FilterIndex]=uigetfile({'*.tif;*.bmp;*.jpg;*.png','Allimagefiles'},'Selecttheinputimage');%获取所选图像的信息（文件名、路径等）ifFilterIndexpic_orig=imread([PathName,FileName]);%读取原始标准图像elsereturnendfigure(1);imshow(pic_orig);%显示原始标准图像title('原始标准图像');[x,y,z]=size(pic_orig);%判断该图像是彩色图像还是灰度图像，如果是彩色图像把它转换为灰度图像if(z~=1)pic_gray=rgb2gray(pic_orig);endfigure(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);%特征提取结果fori=2:x-1forj=2:y-1dg1=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个差分算子进行升序排列iftemp(3)==1%有任意两个大于阈值，则该像素有可能是一特征点result(i,j)=255;%是初选点elseresult(i,j)=0;%不是初选点endendendfigure(4);imshow(result);title('初选点');%第二步：在以初选点为中心的3*3的窗口中计算协方差矩阵与圆度wMatrix=zeros(x,y);%权重矩阵Tq=0.8;%阈值，可设置fori=2:x-1forj=2:y-1ifresult(i,j)==255%如果是初选点gu2=0.0;gv2=0.0;guv=0.0;forii=i-1:iforjj=j-1:jgu2=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));endendDetN=gu2*gv2-guv^2;trN=gu2+gv2;q=4*DetN/(trN*trN);%第三步：设定阈值Tq,若满足则计算权值ifqTqwMatrix(i,j)=DetN/trN;elseresult(i,j)=0;endendendendfigure(5);imshow(result);title('进一步提取的特征点');%第四步：以权值为基础，在一定窗口内抑制局部非最大值候选点,取出局部极大值点vwsize=5;%选择5*5的窗口wradius=floor(vwsize/2);fori=wradius+1:x-wradiusforj=wradius+1:y-wradiusifresult(i,j)==255tempiv=wMatrix(i-wradius:i+wradius,j-wradius:j+wradius);tempsort=sort(tempiv(:),'descend');%将区域内像素按从大至小排列ifwMatrix(i,j)==tempsort(1)&&wMatrix(i,j)~=tempsort(2)%排除整个区域权值相等的情况;elseresult(i,j)=0;endendendendfigure(6);imshow(result);title('forstner提取的最终匹配点');