forstner算子提取特征点文档格式.doc

forstner算子提取特征点文档格式.doc

《forstner算子提取特征点文档格式.doc》由会员分享，可在线阅读，更多相关《forstner算子提取特征点文档格式.doc（4页珍藏版）》请在冰豆网上搜索。

functionForstner,image,vwsize=vwsize,ithresh=ithresh,Tq=Tq

IFN_Elements（vwsize）eq0THENvwsize=5

IFN_Elements（ithresh）eq0THENithresh=50

IFN_Elements（Tq）eq0THENTq=0.5

image=float（image）

imgSize=Size（image,/Dimensions）

xsize=imgSize[0]

ysize=imgSize[1]

灰度的协方差矩阵

result=fltarr（xsize,ysize）

第一步：

利用差分算子提取初选点

fori=1,xsize-2dobegin

forj=1,ysize-2dobegin

dg1=abs（image[i,j]-image[i+1,j]）

dg2=abs（image[i,j]-image[i,j+1]）

dg3=abs（image[i,j]-image[i-1,j]）

dg4=abs（image[i,j]-image[i,j-1]）

dg=[dg1,dg2,dg3,dg4]

temp=dg[sort（dg）]

iftemp[2]gtithreshthenbegin

result[i,j]=255

endifelsebegin

result[i,j]=0

endelse

endfor

endfor

第二步：

在以初选点为中心的3*3的窗口中计算协方差矩阵与圆度

此处可用where提高循环效率

权重矩阵

wMatrix=fltarr（xsize,ysize）

;

是初选点

ifresult[i,j]eq255thenbegin

gu2=0.0&

gv2=0.0&

guv=0.0

forii=-1,1dobegin

forjj=-1,1dobegin

gu2=gu2+（image[i+1,j+1]-image[i,j]）^2

gv2=gv2+（image[i,j+1]-image[i+1,j]）^2

guv=guv+（image[i+1,j+1]-image[i,j]）*（image[i,j+1]-image[i+1,j]）

DetN=gu2*gv2-guv

trN=gu2+gv2

q=4*DetN/（trN*trN）

第三步：

设定阈值Tq,若满足则计算权值

ifqgtTqthenwMatrix[i,j]=DetN/trN

endif

第四步：

以权值为基础，在一定窗口内抑制局部非最大值候选点;

取出局部极大值点

wradius=vwsize/2

fori=wradius,xsize-1-wradiusdobegin

forj=wradius,ysize-1-wradiusdobegin

tempiv=wMatrix[i-wradius:

i+wradius,j-wradius:

j+wradius]

将区域内像素按从大至小排列

tempsort=tempiv（REVERSE（SORT（tempiv）））

排除整个区域像素相等的情况

if（wMatrix[i,j]eqtempsort[0]）and（wMatrix[i,j]netempsort[1]）thenbegin

return,result

end

--------------------

proForstner_test

DEVICE,DECOMPOSED=1

获取本程序所在文件路径

RootDir=Sourceroot（）

file=RootDir+'

\small.bmp'

file=RootDir+'

\8bit_house.bmp'

queryStatus=QUERY_IMAGE（file,imgInfo）

ifqueryStatuseq0thenbegin

Result=DIALOG_MESSAGE（'

参考图像格式不可识别！

'

/error,title='

警告'

）

return

endif

if（imgInfo.CHANNELSne1）thenbegin

图像格式必须为8bit'

imgSize=imgInfo.dimensions

xsize=imgsize[0]

ysize=imgsize[1]

image=READ_IMAGE（file）

resultimg=Forstner（image,ithresh=70,Tq=0.5）

temp=image

index=where（resultimgeq255,count）

print,count

dims=size（resultimg,/dimensions）

ncol=dims[0]

col_index=indexmodncol;

列数

row_index=index/ncol;

行数

WINDOW,/free,XSIZE=xsize*2,YSIZE=ysize

tv,image,0

tv,image,1

fori=0,count-1dobegin

PLOTS,[col_index[i]-2,col_index[i]+2],[row_index[i],row_index[i]],/DEVICE,color='

0000ff'

xl

PLOTS,[col_index[i],col_index[i]],[row_index[i]-2,row_index[i]+2],/DEVICE,color='