matlab编程实现Forstner算子

clear all;
clc;
%获取文件名、路径
[Filename,Pathname,FilterIndex]=uigetfile({'test1.png','D:\新建文件夹\遥感\数字摄影测量\Forstner算子'},'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('灰度图像');
%lever=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);%对四个差分算子进行升序排列
        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提取的最终匹配点');

结果图：