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='警告')
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(resultimg eq 255,count)
print,count
dims=size(resultimg,/dimensions)
ncol=dims[0]
col_index=index mod ncol ;列数
row_index=index/ncol ;行数
WINDOW, /free, XSIZE = xsize*2, YSIZE = ysize
tv,image,0
tv,image,1
for i=0,count-1 do begin
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。