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史振偉（攝影測量作業）\