forstner算子提取特征点

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

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

;Forstner算子

;; image:输入原始图像

; vwsize:窗口宽度

; ithresh:初选差分阈值

; qthresh:兴趣值阈值

function Forstner,image,vwsize=vwsize,ithresh=ithresh,Tq=Tq

IF N_Elements(vwsize) eq 0 THEN vwsize=5

IF N_Elements(ithresh) eq 0 THEN ithresh=50

IF 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

endfor

endfor

;第二步：在以初选点为中心的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

endfor

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

wradius=vwsize/2

for 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

endfor

endfor

return,result

end

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

pro Forstner_test

DEVICE, DECOMPOSED=1

; 获取本程序所在文件路径

RootDir = Sourceroot()

; file=RootDir+'\small.bmp'

file=RootDir+'\8bit_house.bmp'

queryStatus = QUERY_IMAGE(file, imgInfo)

if queryStatus eq 0 then begin

Result = DIALOG_MESSAGE('参考图像格式不可识别！',/error,title='警告')

return

endif

if (imgInfo.CHANNELS ne 1) then begin

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