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进行单波段的输出再进行清晰度的评价，而是通过循环进行每一个波段的遍历。
