Moravec算子和Forstner算子提取图像特征点

采用matlab分别编写的Moravec和Forstner算子提取图像中的特征点，效果不错。

clc;
clear;
tic;
origpic=imread（‘F:\学习资料\N研究生用\研究生课程作业\现代摄影测量\寒假作业moravec&fostner算子\徽.jpg‘）;
origpicgray=rgb2gray（origpic）;    %RGB图像转灰度图像
%%  计算图像长宽像素值mn。
[mn]=size（origpicgray）;   %宽m，长n
%%  利用差分算子（像素（cr）在上下左右四个方向的灰度差分绝对值）提取初选点。
g=double（origpicgray）;
T=30;               %阈值
Initpointcoord=[];   %初选点坐标
for i=2:m-1
    for j=2:n-1
        dg1=abs（g（ij）-g（i-1j））;   %上侧
        dg2=abs（g（ij）-g（i+1j））;   %下侧
        dg3=abs（g（ij）-g（ij-1））;   %左侧
        dg4=abs（g（ij）-g（ij+1））;   %右侧
        dg=[dg1 dg2 dg3 dg4];
        temp=sort（dg）;
        if temp（3）>T
            Initpointcoord=[Initpointcoord;i j];
        end
    end
end
figure（1）;
title（‘初选点‘）
imshow（origpicgray）;
hold on;
plot（Initpointcoord（:2）Initpointcoord（:1）‘r+‘）;
%%  以初选点（cr）为中心的w*w窗口中，计算协方差矩阵N与误差椭圆的圆度q确定备选点。
w=3;                          %窗口宽的值
Initlen=length（Initpointcoord）;   %初选点个数
k=floor（w/2）;
Tq=0.32;                       %阈值Tq=0.32-0.5
candidatepoint=[];             %候选点
for no=1:Initlen
    Initpoint=Initpointcoord（no:）;
    c=Initpoint（1）;
    r=Initpoint（2）;
    h=c-k;
    l=c+k-1;
    p=r-k;
    w=r+k-1;
    gu2=0;gv2=0;guv=0;
    for e=h:l
        for f=p:w
            gu2=gu2+（g（e+1f+1）-g（ef））^2;
            gv2=gv2+（g（ef+1）-g（e+1f））^2;
            guv=guv+（g（e+1f+1）-g（ef））*（g（ef+1）-g（e+1f））;
        end
    end
    DetN=gu2*gv2-guv^2;
    trN=gu2+gv2;
    q=4*DetN/（trN*trN）;
    if q>Tq
        candidatepoint=[candidatepoint;c r];
    end
end
figure（2）;
title（‘备选点‘）
imshow（origpicgray）;
hold on;
plot（candidatepoint（:2）candidatepoint（:1）‘r+‘）;
%%  以权值为依据，选取一适当窗口（w2*w2）中的权值最大者点为特征点。
clear h p q j c r
w2=5;
cadipotlen=length（candidatepoint）;
k2=floor（5/w2）;
featurepoint=[];
h=floor（m/w2）;   %窗口行数
d=floor（n/w2）;   %窗口列数
for p=1:h
    for q=1:d
        tp=（p-1）*w2;
        tq=（q-1）*w2;
        window=zeros（cadipotlen2）;
        for j=1:cadipotlen
            cadipoint=candidatepoint（j:）;
            c=cadipoint（1）;
            r=cadipoint（2）;
           if （c>tp&&ctq&&r               window（j:）=[cr];
           end
        end
        if max（window）~=0
            num=find（window（:1）~=0）;
            pointcoord=window（num:）;
            [pointcoordlen null]=size（pointcoord）;
            weight=[];
            for ii=1:pointcoordlen
                point=pointcoord（ii:）;
                cc=point（1）;
                rr=point（2）;
                h=cc-k;
                l=cc+k-1;
                pp=rr-k;
                w=rr+k-1;
                gu2=0;gv2=0;guv=0;
                for e=h:l
                    for f=pp:w
                        gu2=gu2+（g（e+1f+1）-g（ef））^2;
                        gv2=gv2+（g（ef+1）-g（e+1f））^2;
                        guv=guv+（g（e+1f+1）-g（ef））*（g（ef+1）-g（e+1f））;
                    end
                end
                DetN=gu2*gv2-guv^2;
                trN=gu2+gv2;
         

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        3588  2014-03-02 10:55  2112130006史振伟（摄影测量作业）\Forstner_szw.m
     文件      655569  2014-03-02 11:11  2112130006史振伟（摄影测量作业）\Moravec&Forstner特征提取结果比较.docx
     文件        1752  2014-03-02 10:53  2112130006史振伟（摄影测量作业）\Moravec_szw.m
     文件       20187  2014-02-21 15:17  2112130006史振伟（摄影测量作业）\徽.jpg
     目录           0  2014-03-03 15:27  2112130006史振伟（摄影测量作业）\