01 目的
1. 掌握清晰度的概念;
2. 熟悉清晰度评价指标;
3. 基本掌握清晰度定量指标的 IDL 编程
02 清晰度算法原理介绍
最简单的Brenner函数(梯度滤波器):只需要计算X方向上相差两个像素点的差分,即计算二阶梯度,其公式如下:
接着是能量梯度函数,分别采用三种邻域范围进行图像清晰度评价,公式如下:
(注意: Dist计算时,将单个像元的长宽均视为一个单位计算,如有需要,可读取影像的分辨率进行)
邻域范围如下:
03 程序说明
;基于Brenner函数的图像清晰度评价
;参考文献:冯精武, 喻擎苍, 芦宁, et al. 调焦系统中数字图像清晰度评价函数的研究[J]. 机电工程, 2011, 28(3):354-356.
function way0, data, nl, ns
; 使用Brenner函数计算图像清晰度(邻域为中心像元的中下距离2个像元位置的像元)
DenValue1 = (nl - 1)*(ns - 1)
Tidudata = 0.0
FOR i = 0,ns-3 DO BEGIN
FOR j = 0,nl-3 DO BEGIN
;计算像元梯度值
Tidudata = Tidudata + abs(data[i,j+2]-data[i,j])
ENDFOR
ENDFOR
Tidudata = Tidudata * 1.0 / DenValue1
return, Tidudata
end
function way1, data, nl, ns
; 使用能量梯度函数进行清晰度的计算(邻域为中下、中右)
;定义清晰度计算公式中的分母denominator
DenValue = (nl - 1)*(ns - 1)
Tidudata = 0.0
FOR i = 0,ns-2 DO BEGIN
FOR j = 0,nl-2 DO BEGIN
;计算像元梯度值
Tidudata_1 = abs(data[i,j+1]-data[i,j])
Tidudata_2 = abs(data[i+1,j]-data[i,j])
Tidudata = Tidudata + (Tidudata_1 + Tidudata_2)
ENDFOR
ENDFOR
Tidudata = Tidudata * 1.0 / DenValue
return, Tidudata
end
function way2, data, nl, ns
; 使用能量梯度函数进行清晰度的计算(邻域为十字方向的四个临近像元)
;定义清晰度计算公式中的分母denominator
DenValue = (nl - 2)*(ns - 2)
Tidudata = 0.0
FOR i = 1,ns-2 DO BEGIN
FOR j = 1,nl-2 DO BEGIN
;计算像元梯度值
Tidudata_1 = abs(data[i,j+1]-data[i,j])
Tidudata_2 = abs(data[i+1,j]-data[i,j])
Tidudata_3 = abs(data[i-1,j]-data[i,j])
Tidudata_4 = abs(data[i,j-1]-data[i,j])
Tidudata = Tidudata + (Tidudata_1 + Tidudata_2 + Tidudata_3 + Tidudata_4)
ENDFOR
ENDFOR
Tidudata = Tidudata * 1.0 / DenValue
return, Tidudata
end
function way3, data, nl, ns
; 使用能量梯度函数进行清晰度的计算(邻域为周围八个像元)
;定义清晰度计算公式中的分母denominator
DenValue = (nl - 2)*(ns - 2)
Tidudata=0.0
FOR i = 1,ns-2 DO BEGIN
FOR j = 1,nl-2 DO BEGIN
;计算像元梯度值(也可用循环)
Tidudata_1 = abs(data[i,j+1]-data[i,j])
Tidudata_2 = abs(data[i+1,j]-data[i,j])
Tidudata_3 = abs(data[i-1,j]-data[i,j])
Tidudata_4 = abs(data[i,j-1]-data[i,j])
Tidudata_5 = abs(data[i-1,j-1]-data[i,j]) / sqrt(2.0)
Tidudata_6 = abs(data[i+1,j-1]-data[i,j]) / sqrt(2.0)
Tidudata_7 = abs(data[i-1,j+1]-data[i,j]) / sqrt(2.0)
Tidudata_8 = abs(data[i+1,j+1]-data[i,j]) / sqrt(2.0)
Tidudata = Tidudata + (Tidudata_1 + Tidudata_2 + Tidudata_3 + Tidudata_4 + Tidudata_5 + Tidudata_6 $
+ Tidudata_7 + Tidudata_8)
ENDFOR
ENDFOR
Tidudata = Tidudata * 1.0 / DenValue
return, Tidudata
end
PRO Definition2D
; 01 准备工作
;避免编译时ENVI函数找不到的情形;
COMPILE_OPT IDL2
;初始化ENVI
ENVI,/restore_base_save_files
ENVI_BATCH_INIT
;给出输入文件地址
inputfile = 'D:\task\RemoteSencingImageProcessing_HuiChen\IDL\experiment_3\can_tmr.img'
;将计算的梯度值保存到txt文件中
;定义一个输出文件
outputfile = inputfile+'_Definition2D'+'.txt'
PRINT,outputfile
;OPenw函数为打开一个文件,读出或写入
OPENW,lun,outputfile,/get_lun
; 01 end
; 02 获取数据
;打开Inputfile,返回fid,可以通过该值获得数据的任何信息
ENVI_OPEN_FILE, inputfile, r_fid=fid
;需要判断fid,如果fid返回值为-1,则数据不存在
IF (fid EQ -1) THEN BEGIN
ENVI_BATCH_EXIT
RETURN
ENDIF
;基于fid,利用ENVI_FILE_QUERY获得关于数据文件的信息:包括行列数,波段数,空间维数,文件名,数据类型、数据存储方式等
ENVI_FILE_QUERY, fid, dims=dims, nb=nb ,nl=nl,ns=ns
; 02 end
; 03 计算图像清晰度
;基于ENVI_FILE_QUERY返回的信息,利用ENVI_GET_DATA工具获取数据
;需要注意该工具每次只能获取一个波段, 所以循环获取所有波段并进行处理
for b_i = 0, nb - 1 do begin
data = ENVI_GET_DATA(DIMS=dims, FID=fid , POS=b_i)
data = float(data) ; 将字节类型转化为浮点型(避免后续出现字节与字节型相加出现超出范围导致计算错误的问题)
;定义一个新数组,存储每一个像元的梯度值
;MAKE_ARRAY函数使您能够动态创建一个直到运行时才知道其特征的数组
; Tidudata = make_array(ns,nl-2)
;指定单个像元距离值,此处假设为1
; PixelDist =1.0
; 计算图像清晰度
Tidudata0 = way0(data, nl, ns)
Tidudata1 = way1(data, nl, ns)
Tidudata2 = way2(data, nl, ns)
Tidudata3 = way3(data, nl, ns)
; 输出到命令行窗口
PRINT,Tidudata0
PRINT,Tidudata1
PRINT,Tidudata2
PRINT,Tidudata3
;写入打开的文件
PRINTF,lun, '波段', strcompress(string(b_i)), '的遥感影像Brenner函数清晰度值为',Tidudata0
PRINTF,lun, '波段', strcompress(string(b_i)),'遥感影像能量函数清晰度值为(邻域为2)',Tidudata1
PRINTF,lun, '波段', strcompress(string(b_i)),'遥感影像能量函数清晰度值为(邻域为4)',Tidudata2
PRINTF,lun, '波段', strcompress(string(b_i)),'遥感影像能量函数清晰度值为(邻域为8)',Tidudata3
printf, lun, string(10B) ; 换行 , ASCII中10表示换行符
endfor
; 关闭文件释放资源
ENVI_FILE_MNG, ID = fid, /REMOVE
;释放LUN
FREE_LUN,lun
; 03 end
; 当用户使用ENVI批处理功能批量处理大量数据时,可以使用ENVI_BATCH_EXIT函数来关闭所有打开的文件、删除所有临时文件并释放内存。
ENVI_BATCH_EXIT ; 反正都程序结束了,用不用都无所谓
END
04 总结与感想
这次实验,主要有几个难点,一个是难以想象需要将矩阵data转化为float(因为无法想象这会造成后续计算的范围超出导致计算错误);第二个编写将各个计算图像清晰度的方法,但是实际上并不难,只是稍显繁琐;第三个就是对每一个波段影像对进行清晰度的评价,这里并不使用ENVI进行save as进行单波段的输出再进行清晰度的评价,而是通过循环进行每一个波段的遍历。