forstner算子提取特征点.doc

Forstner算子提取特征点（原创）;-;Forstner算子; image:输入原始图像; vwsize:窗口宽度; ithresh:初选差分阈值; qthresh:兴趣值阈值function Forstner,image,vwsize=vwsize,ithresh=ithresh,Tq=TqIF N_Elements(vwsize) eq 0 THEN vwsize=5IF N_Elements(ithresh) eq 0 THEN ithresh=50IF N_Elements(Tq) eq 0 THEN Tq=0.5image=float(image)imgSize = Size(image, /Dimensions)xsize=imgSize0ysize=imgSize1;灰度的协方差矩阵result=fltarr(xsize,ysize);第一步：利用差分算子提取初选点for i=1,xsize-2 do begin for j=1,ysize-2 do begin dg1=abs(imagei,j-imagei+1,j) dg2=abs(imagei,j-imagei,j+1) dg3=abs(imagei,j-imagei-1,j) dg4=abs(imagei,j-imagei,j-1) dg=dg1,dg2,dg3,dg4 temp=dgsort(dg) if temp2 gt ithresh then begin resulti,j=255 endif else begin resulti,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 resulti,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+(imagei+1,j+1-imagei,j)2 gv2=gv2+(imagei,j+1-imagei+1,j)2 guv=guv+(imagei+1,j+1-imagei,j)*(imagei,j+1-imagei+1,j) endfor endfor DetN=gu2*gv2-guv trN=gu2+gv2 q=4*DetN/(trN*trN) ;第三步：设定阈值Tq,若满足则计算权值 if q gt Tq then wMatrixi,j=DetN/trN endif endforendfor;第四步：以权值为基础，在一定窗口内抑制局部非最大值候选点;取出局部极大值点wradius=vwsize/2for i=wradius,xsize-1-wradius do begin for j=wradius,ysize-1-wradius do begin tempiv=wMatrixi-wradius:i+wradius,j-wradius:j+wradius ;将区域内像素按从大至小排列 tempsort=tempiv(REVERSE(SORT(tempiv) ;排除整个区域像素相等的情况 if (wMatrixi,j eq tempsort0) and (wMatrixi,j ne tempsort1) then begin resulti,j=255 endif else begin resulti,j=0 endelse endforendforreturn,resultend;-pro Forstner_testDEVICE, DECOMPOSED=1; 获取本程序所在文件路径RootDir = Sourceroot(); file=RootDir+small.bmpfile=RootDir+8bit_house.bmpqueryStatus = QUERY_IMAGE(file, imgInfo)if queryStatus eq 0 then begin Result = DIALOG_MESSAGE(参考图像格式不可识别！,/error,title=警告) returnendifif (imgInfo.CHANNELS ne 1) then begin Result = DIALOG_MESSAGE(图像格式必须为8bit,/error,title=警告) returnendifimgSize = imgInfo.dimensionsxsize=imgsize0ysize=imgsize1image=READ_IMAGE(file)resultimg=Forstner(image,ithresh=70,Tq=0.5)temp=imageindex=where(resultimg eq 255,count)print,countdims=size(resultimg,/dimensions)ncol=dims0col_index=index mod ncol ;列数row_index=index/ncol ;行数WINDOW, /free, XSIZE = xsize*2, YSIZE = ysizetv,image,0tv,image,1for i=0,count-1 do begin PLOTS,col_indexi-2,col_indexi+2,row_ind