forstner算子提取特征点.docVIP

Forstner算子提取特征点（原创） ;------------------------------;Forstner算子 ; 参考自：张祖勋，数字摄影测量学。P125; image:输入原始图像; vwsize:窗口宽度; ithresh:初选差分阈值; qthresh:兴趣值阈值function Forstner,image,vwsize=vwsize,ithresh=ithresh,Tq=Tq IF N_Elements(vwsize) eq 0 THEN vwsize=5IF N_Elements(ithresh) eq 0 THEN ithresh=50IF N_Elements(Tq) eq 0 THEN Tq=0.5 image=float(image) imgSize = Size(image, /Dimensions)xsize=imgSize[0]ysize=imgSize[1] ;灰度的协方差矩阵result=fltarr(xsize,ysize) ;第一步：利用差分算子提取初选点 for i=1,xsize-2 do begin?? for j=1,ysize-2 do begin??? 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)] ??? if temp[2] gt ithresh then begin???? result[i,j]=255??? endif else begin???? result[i,j]=0??? endelse ?? endforendfor ;第二步：在以初选点为中心的3*3的窗口中计算协方差矩阵与圆度;此处可用where提高循环效率 ;权重矩阵wMatrix=fltarr(xsize,ysize) for i=1,xsize-2 do begin?? for j=1,ysize-2 do begin??? ;是初选点??? if result[i,j] eq 255 then begin???? gu2=0.0 gv2=0.0 guv=0.0???? for ii=-1,1 do begin????? for jj=-1,1 do begin?????? 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,若满足则计算权值???? if q gt Tq then wMatrix[i,j]=DetN/trN ??? endif ?? endforendfor ;第四步：以权值为基础，在一定窗口内抑制局部非最大值候选点 ;取出局部极大值点wradius=vwsize/2for i=wradius,xsize-1-wradius do begin?? for j=wradius,ysize-1-wradius do begin??? tempiv=wMatrix[i-wradius:i+wradius,j-wradius:j+wradius]??? ;将区域内像素按从大至小排列??? tempsort=tempiv(REVERSE(SORT(tempiv)))??? ;排除整个区域像素相等的情况??? if (wMatrix[i,j] eq tempsort[0]) and (wMatrix[i,j] ne tempsort[1]) then begin???? result[i,j]=255??? endif else begin???? result[i,j]=0??? endelse?? endforendfor return,resultend ;--------------------pro Forstner_testDEVICE, DECOMPOSED=1 ; 获取本程序所在