forstner算子提取特征点.doc

forstner算子提取特征点.doc

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

Forstner算子提取特征点（原创）

;------------------------------

;Forstner算子

;;image:

输入原始图像

;vwsize:

窗口宽度

;ithresh:

初选差分阈值

;qthresh:

兴趣值阈值

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）

fori=1,xsize-2dobegin

  forj=1,ysize-2dobegin

   ;是初选点

   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]）

     endfor

    endfor

    DetN=gu2*gv2-guv

    trN=gu2+gv2

    q=4*DetN/（trN*trN）

    ;第三步：

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

    ifqgtTqthenwMatrix[i,j]=DetN/trN

   endif

  endfor

endfor

;第四步：

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

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

    result[i,j]=255

   endifelsebegin

    result[i,j]=0

   endelse

  endfor

endfor

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

  Result=DIALOG_MESSAGE（'图像格式必须为8bit',/error,title='警告'）

  return

endif

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='0000ff'xl

endfor

end