大文件遥感图像的中值滤波算法研究.doc

大文件遥感图像的中值滤波算法研究唐山师范学院计算机与科学系 09级本科班 翟艳侠 091220141081摘要: 利用中值滤波技术能够很好的抑制噪声和保护边缘信息这一特性, 用不同的滤波窗口对遥感图像进行处理, 并做出处理后图像的归一化差值图像。通过图形分析发现差值图像上的信息比原始图像细节信息更丰富,但是当遇到比较大的文件，我们就需要对中值滤波算法进行优化，使它能更好的为我们服务。关键词: 中值滤波; 图像处理; 遥感; 遥感图像处理Abstract: the use of median filtering technology can very good suppression noise and protection of the edge information of the properties, use different filtering window in remote sensing image processing, image processing and make after the normalized difference image. Through the graphic analysis found that difference image information than the original image detail information more rich, but when meet larger files, we need to median filtering algorithm for optimization, and make it better for our service.Keywords: median filtering; Image processing; Remote sensing; Remote sensing image processing一、引言遥感图像作为遥感技术的重要成果, 对各个领域的决策及现象分析起到了很大的作用 。 所以在使用一幅遥感图像之前, 必须对遥感图像进行处理 。目前, 常用的图像噪声滤波方法有低通滤波法 、带通滤波法 、 平滑 、 锐化滤波 、 均值滤波等, 但是这些方法在消除噪声的同时也会将有用的信息一并消除, 影响了图像的真实性, 所以很多人试图用中值滤波法对图像进行处理 。 研究表明, 中值滤波法是消除随机图像噪声的最佳方法, 对图像进行中值滤波后,可以很好地消除分布在图像上的各类随机噪声 。中值滤波 ( M F ) 是一项比较成熟的技术 。 早在20 世纪70 年代, L. R. R ab m e r, M . R. Sa m b u r 和C . E . Sch m i d t 就提出了早期的中值滤波算法原形 。 在1982 年, T. A. N o de s 和 N. C . Ga l l agh e r 对中值滤波技术进行了分析并对其特性进行了描述, 为以后中值滤波的发展奠定了基础 。 1986 年 G . R. A rce, N. C . Ga l l agh e r 和T. A. N o ed s 进一步总结了它的理论基础和应用领域, 使中值滤波在信号处理的领域得到重视 。 国内较早对中值滤波技术进行研究的是1985 年由李承楚发表的 “关于中值滤波的理论基础” 一文 。 虽然国内中值滤波技术研究起步比较晚,但是对于新方法的研究却发展较快 。由于中值滤波同时具有抑制噪声 ( 尤其是脉冲噪声) 和保护边缘的特征, 而在信号处理领域得到很大的重视, 特别在非平稳信号的处理中取得了较大的成功, 随之而来, 各种中值滤波的改进技术也在不断的发展, 应用领域在不断的拓宽 。二、中值滤波算法过程本文提出的中值滤波的快速算法的基本思想是：原始数据序列上中值滤波的滑窗在移动过程中，当前窗只要删除其最早的元素，加入窗后的新元素，即成为下一窗的内容。下一窗的中值滤波实现可利用上次中值滤波的排序结果，新元素的插人位置用有序序列快速查找算法求得，新元素插人与最早的元素删除的实现采用独特的数据结构，将是新元素覆盖最早的元素，即是插人兼并了删除。设置(2N+1)个连续存贮单元(存放浮点数) (。)t (1)， (2)组成的循环序列用来存放窗内的数据元素 按照先进先出的原则，后来的数据元素总是替换当前最早存放的数据元素。设置(2+1)个连续存贮单元(存放整数) (。)， (1)，s(2N)顺序存放的是，若上述窗内元素从小至大排序后，顺序的元素在 序列中的下标值，即满足 ( (。)w(s(1) w( (2)。设置(2+1)个连续存贮单元(存放整数)a(0)，n(1)，a(ZN)分别存放 s 序列中存有其下标的存贮单元的下标值，即满足 (。(f)， =O，1，2N。可以这样认为，把 序列和n 序列中具有同一下标的两个存贮单元当作独立结构单元，s 序列中一存贮单元指向上述某一结构单元，这个结构单元中的。存贮单元值表示了这个结构单元指向该 s 存贮单元。下面实现中值滤波的快速算法。首先，令 序列中的存贮单元值全为零，s 序列和 n 序列中的存贮单元分别存放各自的下标值，即 (f)=0， (f)=f，n(f)=f，f=O，1 。,2N。另外，设置下标 =0第一步，求当前准备进入窗的数据元素 x(ra)在 s 序列中的内插位置，用对半查找算法实现脚，如图 1 所示。在图 1 中，有序序列对半查找的区问下界为工，上界为 h，中部为 ，通过比较 (s(1)与待查量 (m)，若不相等，则调整 L 或 h，使下次查找的区问比前次的减少一半。输出的 反映了 (s(1) (m) (s( +1)，其中：一 1 2N。由于 (s(一 1)和(s(2+1)不存在，可把它们当作特殊的符号，假想 (s(一 1)为无穷小， (#(2+1)为无穷大。对由 n 个元素组成的有序序列，若是用对半查找算法，至多做与 l+1og：n个元素的比较就能找到元素的内插位置 第二步，在 W，s，n 序列中插入 x(ra)的有关信息，如图 2 所示。若n( )，则将 s 序列第l+1 至第 n( )一 1 问各存贮单元所存值顺序地移到下一单元，并将该单元所存值作为下标对应的。序列中存贮单元指向该 s 单元，最后赋值 ( )= (m)，n( )=l+1，s(k+1)= 若。( ) ，则将 s序列中第 n()+1 至第 问各存贮单元所存值顺序地移到前一单元，并要将该单元所存值作为下标对应的 n 序列中存贮单元指向该 s 单元，最后赋值 ( )( )，口( )=z,5(1)= 。此时，w(s()为序列x(m 一 2N)， ( 一 2N+1)， (m)的中值。第三步，令，为( +1)rood(2+1)，m 增加 1，取下一个 值，返回第一步按同样步骤求下一中值。如此反复进行，快速求得一系列中值。S 和 n 序列中存贮单元的大小若是选择为单字节(计算机内存最小分配单元)，可保证中值滤波窗口长度达到 2s；2s6，这能满足大多数中值滤波的要求。于是，图 2 中的针对存贮单元的操作显得快捷。总之，本文提出的中值滤波算法的快速特性在于 t 为了求一个中值，运算时间主要用在做一次对半查找，至多与 1+1ogz(2N+1)个数据元素比较。由于数据元素为任意大小或精度的浮点数，所以窗内任意两个元素完全相同的可能性极小，求元素内插位置时数据元素的比较个数基本上为 1+r-logt(2+1)。另外，运算时间还包括 S 和 n 序列完成数据元素内插时所进行的线性过程。三、遥感图像文件1.读取BSQ的源代码：clear allclclines=400;samples=640;N=6;img=fopen(D:sample_BSQ,rb);for i=1:N bi=fread(img,lines*samples,uint8); band_cov=reshape(bi,samples,lines); band_cov2=band_cov; band_uint8=uint8(band_cov2); tif=imadjust(band_uint8); mkdir(D:MATLAB,tifbands1) name=D:MATLABtifbands1tif,int2str(i),.tif; imwrite(tif,name,tif); tilt=波段,int2str(i); subplot(3,2,i),imshow(tif);title(tilt); endfclose(img);2.读取BIP源代码clear allclclines=400;samples=640;N=6;for i=1:N img=fopen(D:MATLABsample_BIP,rb); b0=fread(img,i-1,uint8); b=fread(img,lines*samples,uint8,(N-1); band_cov=reshape(b,samples,lines);band_cov2=band_cov;% band_uint8=uint8(band_cov2); tif=imadjust(band_uint8); mkdir(E:MATLAB,tifbands) name=E:MATLABtifbandstif,int2str(i),.tif; imwrite(tif,name,tif); %imwrite(A,filename,fmt) tilt=波段,int2str(i); subplot(3,2,i),imshow(tif);title(tilt); fclose(img);end3.读取BIL的源代码clear allclc lines=400;samples=640;N=6;for i=1:N bi=zeros(lines,samples); for j=1:samples img=fopen(D:MATLABsample_BIL,rb); bb=fread(img,(i-1)*640,uint8); b0=fread(img,1*(j-1),uint8); bandi_linej=fread(img,lines,uint8,1*(N*samples-1); fclose(img); bi(:,j)=bandi_linej; end band_uint8=uint8(bi); tif=imadjust(band_uint8); mkdir(D:MATLAB,tifbands) name=D:MATLABtifbandstif,int2str(i),.tif; imwrite(tif,name,tif); tilt=,int2str(i); subplot(3,2,i),imshow(tif);title(tilt);end4、运行结果 四、中值滤波算法优化1 )在噪声滤波之前 , 先将滤波范围缩小.根据含噪图像的灰度直方图 , 将灰度值按区域分开 , 设定一个阈值 , 把 255 - , 255 与 0, 作为下一步要滤波的范围 , 将含噪声点的区域与信号分开 , 并 且 只 对 含 噪 声 点 的 区 域进行滤波处理. 这样可以有效避免在中值滤波时对所有像素均用滤窗内的中值来代替所造成的图像模糊 , 同时又能节省运算时间.2 )对选定的滤波范围进行噪声检测.设滤窗内所有像素点构成的集合为Sxy= f ( i + s , j + t) | ( i , j ) ( ( 0,) ( 255 - , 255 ) ) , s , t ( 1, 3 ) 其中 , Zm in是 Sxy滤窗内灰度最小值 ; Zmax 是 Sxy滤窗内灰度最大值 ;favg 是 Sxy滤窗内除去最大值 、 最小值后的灰度平均值. 另设定一个检测阈值 T, 如果 f ( i , j )与 favg 的差值大于 T, 而且 f ( i , j ) 与 Zmax, Zm in, ( Zmax- 1 )或 ( Zm in+ 1 )中的某一个值相等 , 则判定 f ( i , j )为噪声点 , 即噪声点满足的条件为 ( f ( i , j ) = Zmax f ( i , j ) = ( Zmax - 1 ) f ( i , j ) =Zm in f ( i , j ) = ( Zm in+ 1 ) ) I ( | f ( i , j ) - favg| E T )其中 , T值对噪声检测的效果有较大影响 , 其值应根据图像中噪声的干扰程度而定. 实验发现 , 噪声干扰小时 , T取较大值 , 反之 , T取较小值.设 f ( x, y )是点 ( x, y )上的灰度值 , Sxy是点 ( x, y )所对应的掩模窗口 , 具体算法实现如下 :1 ) 结合灰度直方图取阈值 和 T ;2 )如果 f ( x, y ) 0且 A2 0 , 则跳转至 5 ) ;4 )增加滤窗 Sxy的尺寸 , 如果 Sxy F Smax, 则跳转至 3 ) , 否则 , 跳转至 6 ) ;5 ) B1= f ( x, y ) - Zm in - 1, B2= Zmax - f ( x, y ) - 1,如果 B1 0 且 B2 0, 则输出 f ( x, y ) , 否则 , 继续执行;6 )如果 | f ( i , j ) - favg| T, 输出 f ( x, y ) ; 否则 ,输出 f ( x, y ) .该算法 在 滤 波 之 前 先 判 断 像 素 点 是 否 在 滤波范围内 , 这 样 可 以 节 省 滤 波 时 间 , 也 可 以 避 免大面积滤 波 造 成 的 图 像 模 糊. 对 于 属 于 滤 波 范围内的像素点 , 须进一步判断滤窗内的中值是否为脉冲点 , 若是脉冲点 , 则增大滤波窗口重新判断 , 只要没有达 到 滤 窗 的 尺 寸 极 限 , 就 一 直 循 环 判 断 ,直到找出中值 ; 若达到了所设定的滤窗极限仍然没有找到中值 , 则通过判断该点的灰度值 f ( x, y ) 与均值 favg 的差是否达到阈值 T 来决定是否输出. 若找到了中值 , 则判断该像素点是否满足噪声点条件 :若满足 , 则输出中值 ; 若不满足 , 则输出该像素点.五、总结本文讨论了标准中值滤波在不同滤波窗口时的滤波情况 ,