朱正为8数字图像处理的应用.docx

1、朱正为8数字图像处理的应用第八章 数字图像处理的应用 本章通过几个MATLAB数字图像处理实例，介绍数字图像处理技术的应用，具体内容包括： MATLAB图像处理工具箱和基本的图像处理操作 图像直方图及其拟合 图像平滑：中值滤波去噪 图像复原8.1 MATLAB图像处理工具箱和基本的图像处理操作图像处理工具箱是一个函数的集合，它扩展了matlab数值计算环境的能力。这个工具箱支持了大量图像处理操作，包括：空间图像变换 Spatial image transformations形态操作 Morphological operations 邻域和块操作 Neighborhood and block o

2、perations线性滤波和滤波器设计 Linear filtering and filter design格式变换 Transforms 图像分析和增强 Image analysis and enhancement图像登记 Image registration清晰化处理 Deblurring兴趣区处理 Region of interest operations 工具箱里的函数都是M文件，可以通过type function_name来查看代码，也可以通过写自己的matlab函数来扩展工具箱。要查看是否安装了该工具箱，可以使用ver命令来查看已安装的工具箱。 其他相关工具箱有：DSP Block

3、setImage Acquisition ToolboxMapping ToolboxSignal Processing ToolboxWavelet Toolbox 下面通过简单的例子介绍工具箱的使用方法，此例子介绍了一些基本的图像处理操作，包括读、写图像，演示图像的直方图均衡，得到图像信息等。（1）MATLAB中的数字图像表示在MATLAB中，一幅数字图像表示成矩阵(以二维数组的形式存储的)：其中的每一个元素称为象素，符号f(p,q)表示位于第p行和q列的元素。（2）读取一个图像使用imread可以将图像读入MATLAB环境。clear, close all f=imread(D:myim

4、agesdianluban.jpg);（3）显示一个图像imshow(f) /显示图像f。imshow(f, G) /G是显示图像的灰度级数。省略时默认灰度级数为256。imshow(f, low, high)/将所有小于或等于low的值都显示为黑色，所有大于或等于high的值都显示为白色。imshow(f, )/将变量low设置为数组f的最小值，将变量high设置为数组f的最小值，这一形式在显示一幅动态范围较小的图像或既有正值又有负值的图像时很有用。 imshow(f), figure, imshow(g)保持第一幅图像并同时显示第二幅图像。（4）检查图像在工作区中的表现形式 whosNam

5、e Size BytesClassI 291x240 69840uint8 arrayGrand total is 69840 elements using 69840 bytes图像数据是按数组的形式存储的。（5）显示图像直方图和图像直方图均衡处理显示图像的直方图 figure, imhist(I)或figure, imhist(I,b)/I是图像文件; b是灰度级的个数，省略默认为256执行直方图均衡化I2=histeq(I);figure, imshow(I2) 均衡后的直方图为figure, imhist(I2) （6） 保存图像（写图像到硬盘） imwrite(f, filename

6、, jpg)或imwrite(f, C:pout2.png)（7）得到图像信息imfinfo(C:pout2.png) / 返回关于图像文件的信息8.2 图像直方图及其拟合（1）概述雷达目标和杂波幅度分布特性是雷达目标和杂波的主要统计特性之一，对雷达信号处理、检测、识别、仿真及对雷达系统设计和性能评估均有十分重要的意义。长期以来，雷达工作者一直在研究和探讨这一问题，由于雷达杂波比较复杂，它包括地物杂波、海杂波、气象杂波、鸟群杂波（仙波）和箔条杂波等各种有源和无源干扰，并且在不同条件下，又是千变万化的，故分析起来比较困难。一般都是采用统计的方法，对它们进行分析或对实测数据进行拟合来对雷达杂波的幅

7、度分布进行建模。到目前为止，雷达杂波幅度分布模型有高斯分布、瑞利分布、K分布、韦布尔分布、对数-正态分布、莱斯（Rice）分布模型等。现代目标模型有莱斯模型、对数-正态分布模型和（Chi方）模型等。本节针对高分辨率SAR图像（的直方图），分析几种常用模型的建模能力（其核心是各分布概率密度函数参数的估计，其中常用的参数估计方法有矩估计法、极大似然估计法等）。（2）需拟合的直方图草地SAR图像、森林SAR图像和城市SAR图像。 （3）举例采用对数正态分布对高分辨率SAR图像的直方图进行拟合。（A）对数正态分布的定义定义：如果随机变量的函数服从正态分布，则称服从参数为和的对数正态分布，记作。对数正态

8、分布的概率密度函数为：（B）参数与的极大似然估计根据极大似然估计法的原理，由上式，得似然函数为：两边取对数，得：似然方程组为：故得与的极大似然估计分别是其中是的无偏估计，是的有偏（渐进无偏）估计。可以证明的无偏估计为：（C）MATLAB程序（a）求对数正态分布参数的mu（）和delta2（）clear all;i1=imread(D:image003.jpg);i2=rgb2gray(i1); %将真彩色图像转换成灰度图像x=double(i2); %创建原图像的副本并将其元素转换为double型，uint8型数据无法进行计算p=size(x,1);q=size(x,2);N=p*q;s=0;

9、%求mu for m=1:pfor n=1:qif x(m,n)=0y(m,n)=0;s=s+1;elsey(m,n)= log(x(m,n); endendendmu=sum(sum(y)/N;%求sigma2for i=1:pfor j=1:qz(i,j)=(y(i,j)-mu).2; endendsigma2= sum(sum(z)/N;musigma2mu=4.8947, 4.6667, 3.8266sigma2=0.0095, 0.0341, 0.0677（b）绘制图像灰度直方图clear all;a=imread(D:image002.jpg);p,q=size(a);x=1:1:

10、256; %设定X轴的范围total=zeros(1,256); %用于存放0255各级灰度出现的次数。zeros(m,n)将产生一个mxn的以零为项目之矩阵；而zeros(n)会产生一个nxn的以零为内容之方矩阵b=double(a); %将创建原图像的副本并将其中的数据元素转换为double型，uint8型的数据无法进行计算sum=0;for i=1:p*q %统计0255各级灰度在图像中出现次数total(b(i)+1)=total(b(i)+1)+1; %由于total的下标从1开始到256而b(i)的范围是0255，所以需要将b(i)的值加1才可以使b(i)和total对应end;f

11、or j=1:256 %统计0255各级灰度的累计出现次数 sum=total(j)+sum;jtem(j)=sum;end;%以下是结果显示部分y=total(x)/(p*q);plot(x,y, k-);axis(0,450,0,(max(total)+100)/ (p*q); %格式axis(xmin xmax ymin ymax)（c）绘制对数正态分布概率密度函数clear all;%mu=4.8947;%sigma2=0.0095;mu=4.6667;sigma2=0.0341;%mu=3.8266;%sigma2=0.0677;sigma=sqrt(sigma2);y=lognpd

12、f(1:300,mu,sigma);hold on;plot(y, b-);xlabel(灰度幅值);ylabel(概率密度);8.3 图像平滑：中值滤波去噪（1）在图像中加入椒盐噪声 f=imread(E:pro7.jpg); % The image pro7.jpg is located in the root directory of E. type = salt & pepper; p = 0.2; RN = imnoise(f, type, p); c = find(RN = 1); %将所有数组RN的索引返回c中，它指向非零元素 fs = f; fs(c) =255; %在图像中加

13、入盐噪声 d = find(RN = 0); fsp = fs; fsp(d) = 0; %在含盐图像中加入椒噪声 figure, imshow(fsp)（2）在图像中加入高斯噪声 f=imread(E:pro7.jpg); type = gaussian; u = 0; var = 1; faddedgn = imnoise(f, type, u, var); imshow(f), figure, imshow(faddedgn)（3）椒盐噪声污染图像中值滤波结果 h31 = medfilt2(fsp,3 3,symmetric); imshow(h31) h32 = medfilt2(h3

14、1,3 3,symmetric); figure, imshow(h32) h33 = medfilt2(h32,3 3,symmetric); figure, imshow(h33) (a) 一次3*3 中值滤波结果 (b) 两次3*3 中值滤波结果 (c) 三次3*3 中值滤波结果（4）高斯噪声污染图像中值滤波结果 f=imread(E:pro7.jpg); type = gaussian; u = 0; var = 1; faddedgn = imnoise(f, type, u, var); imshow(f), figure, imshow(faddedgn) hg1 = medfi

15、lt2(faddedgn,7 7,symmetric); figure, imshow(hg1) hg2 = medfilt2(hg1,7 7,symmetric); figure, imshow(hg2) hg3 = medfilt2(hg2,7 7,symmetric); figure, imshow(hg3) (a) 一次7*7 中值滤波结果 (b) 两次7*7 中值滤波结果 (c) 三次7*7 中值滤波结果8.4 图像复原参阅冈萨雷斯数字图像处理（MATLAB版）5.5-5.10节（P.123-134）8.4.1退化函数的建模成像设备实验法用数学方法把PSF模型化典型方法：通过产生PS

16、F及测试各种复原算法的结果来试验采用盲去卷积来推断PSF本节讨论使用函数imfilter和fspecial，以及噪声生成函数imnoise建模PSF的技术。由场景和传感器产生的图像模糊可用低通滤波器来建模；运动模糊使用IPT（图像处理工具箱）函数fspecial（生成线性空间滤波器）建模。（1）产生运动模糊的PSF:PSF=fspecial(motion, len, theta)调用fspecial将返回PSF，它近似于摄像机有len个像素的线性移动的效果。参数theta以度为单位，以顺时针方向相对正水平轴进行度量。（2）创建一个已知PSF的退化图像，使用函数imfilter（线性空间滤波）：

17、 g=imfilter(f, PSF, circular);其中circular为边界选项，用来减少边界效应。（3）使用函数checkerboard产生测试图案C=checkerboard(NP, M, N)其中NP是每个小正方形一边的像素数，M是行数，N是列数。测试板左半部分的亮正方形是白色的，右半部分的亮正方形是灰色的。若省略了N，则默认为M，若都省略，则默认为8。（4）由于有些复原算法对于大图像很慢，因此通常对小图像做实验，而通过像素复制以放大显示图像（自定义函数）：B=pixeldup(A, m, n)该函数将A的每个像素在垂直方向上总共复制了m次，在垂直方向上总共复制了n次。自定义函

18、数pixeldup(A,m,n)的程序function B=pixeldup(A,m,n)%PIXELDUP duplicates pixels of an image in both directions%Check inputsif nargin2 error(At least two inputs are required.);endif nargin=2 n=m;end%Generatr a vector with elements 1:size(A,1).u=1:size(A,1);%Duplicate each element of the vector m times.m=rou

19、nd(m);%Protect against nonintergers.u=u(ones(1,m),:);u=u(:);%Now repeat for the other direction.v=1:size(A,2);%Duplicate each element of the vector m times.n=round(n);%Protect against nonintergers.v=v(ones(1,n),:);v=v(:);B=A(u,v);（5）举例：模糊噪声图像的建模f=checkerboard(8);/产生测试图像(a)PSF=fspecial(motion, 7, 45)

20、;gb=imfilter(f, PSF, circular);/产生退化图像(b)noise=imnoise(zeros(size(f), gaussian, 0 ,0.001);/产生噪声模式(c)，通常我们直接使用imnoise(gb, gaussian, 0 ,0.001)将噪声加到gb上g=gb + noise;/产生模糊噪声图像(d)imshow(pixeldup(f, 8), );/放大到512*512显示8.4.2维纳滤波复原维纳滤波复原是一种线性图像复原方法。我们首先定义平均噪声功率和平均图像功率及它们的比率：有时用来近似代替，以便产生一个常数数组，产生参数维纳滤波器，可以对直

21、接逆滤波产生重大改进。（1）在IPT中，维纳滤波使用函数deconvwnr(g, PSF)来实现：fr=deconvwnr(g, PSF) /假设噪信比为零，这种形式实现的即是逆滤波器。fr=deconvwnr(g, PSF, NSPR) /假设噪信功率比已知，这用于实现参数维纳滤波器。fr=deconvwnr(g, PSF, NACORR, FACORR) /假设噪声和未退化图像的自相关函数NACORR和 FACORR是已知的。由可知，通过计算图像（或噪声）功率谱的傅里叶逆变换可得到图像（或噪声）的自相关函数。（2）举例：使用函数deconvwnr复原前述（运动）模糊噪声图像fr1=deco

22、nvwnr(g, PSF) /直接逆滤波，结果如图(b)所示，这个结果由噪声的效果所决定。 fr2=deconvwnr(g, PSF, R) /使用比率R复原图像，结果如图(c)所示。比率R是利用噪声图像和原图像得到的，过程如下：Sn=abs(fft2(noise).2;nA=sum(Sn (:)/prod(size(noise);Sf=abs(fft2(f).2;fA=sum(Sf (:)/prod(size(f);R=nA/fA;使用自相关函数和deconvwnr进行图像复原：NCORR=fftshift(real(ifft2(Sn);ICORR=fftshift(real(ifft2(S

23、f);Fr3=deconvwnr(g, PSF, NCORR, ICORR);结果如图(d)所示，所得结果虽然仍有一些噪声存在，但已和原图像很接近，因为原图像和噪声都是已知的，所以可以正确地估算参数。在实践中，当这些量之一（或更多）未知时，挑战便是在实验中智能地选择所用的函数，直到获得可接受的结果为止。8.4.3 约束的最小二乘方滤波复原约束最小二乘方滤波复原也是一种线性图像复原方法。（1）约束最小二乘方滤波复原在IPT中通过函数deconvreg来实现：fr=deconvreg(g, PSF, NOISEPOWER, RANGE)/ NOISEPOWER 与成正比，其较好的初始估计为，最后的

24、结果可能有很大的不同。RANGE为值的范围，在求时，该算法受这个范围的限制，默认范围是10-9,109，MATLAB中的符号为1e-10, 1e10。若将上述两个参数排除在参量之外，则函数deconvreg就会产生一个逆滤波方案。（2）举例：使用函数deconvreg复原前述（运动）模糊噪声图像因为图像大小为64*64，噪声的方差为0.001，均值为0，所以NOISEPOWER的初始估计为。图（A）显示了使用下述命令得到的结果：fr=deconvreg(g, PSF, 4); 其中g和PSF来自8.4.1节例子。可以看出NOISEPOWER的值并不是非常好，通过对参数NOISEPOWER和RA

25、NGE进行实验，得出如图（B）所示的结果，它是使用如下命令得到的：fr=deconvreg(g, PSF, 0.4, 1e-7, 1e7);NOISEPOWER的值下调了一个数量级，而RANGE比默认值更加紧缩。维纳滤波结果比该结果要好得多，但其条件是噪声和图像谱已知，若没有这些信息，则用两种滤波器通过实验得到的结果常常是可比的。8.4.4 使用Lucy-Richardson算法的迭代非线性复原（1）L-R算法实现作为复原的工具，非线性迭代技术常常会获得比线性方法更好的结果。在IPT工具箱中，选择的非线性方法是由Richardson和Lucy开发的技术（L-R算法是从最大似然公式引出来的），L

26、-R算法由函数deconvlucy实现：fr=deconvlucy(g, PSF, NUMIT, DAMPAR, WEIGHT)其中参数NUMIT为迭代次数（默认为10次）。参数DAMPAR指定结果图像与原图像g之间的偏离阈值。默认值为0。当像素值偏离原值的范围在DAMPAR之内时，就不用再迭代。WEIGHT是一个与输入图像g同样大小的数组（默认值是与g同样大小的一个单位数组），它为每一个像素分配一个权重来反映其逼真的程度。例如一个不良像素会被赋以权重值零，从而paic 该像素来求解。当用一个指定的PSF来模拟模糊时，WEIGHT可从计算像素中剔除那些来自图像边界的像素点，若PSF的大小为n*

27、n，则在WEIGHT中用到的零边界宽度为ceil(n/2)。（2）举例f=checkerboard(8);/产生测试图像(a)imshow(pixeldup(f, 8); /放大显示PSF=fspecial(gaussian, 7, 10);/产生一个大小为7*7且标准偏差为10的高斯（低通滤波模糊）PSFSD=0.01;g=imnoise(imfilter(f, PSF), gaussian, 0, SD.2);/用PSF产生高斯低通模糊的退化图像，同时添加均值为0、标准偏差为0.01的高斯噪声/余下部分使用函数deconvlucy对图像g进行复原处理DAMPAR=10*SD; /将DAMP

28、AR赋以10倍于SD的值；LIM=ceil(size(PSF, 1)/2);/产生数组WEIGHT，WEIGHT数组的大小是64*64，且有值为0的4像素宽的边界，其余的像素都是1WEIGHT=zeros(size(g);WEIGHT(LIM+1:end- LIM, LIM+1:end- LIM)=1;NUMIT=5；f5=deconvlucy(g, PSF, NUMIT, DAMPAR, WEIGHT);imshow(pixeldup(f5, 8); /放大显示(a)原图像；(b)由高斯噪声污染和高斯低通模糊的图像；(c)到(f)为对图像(b)用L-R算法分别迭代5次、10次、20次和100

29、次后的复原图像。在所有结果图像中，细黑边界都是由数组WEIGHT中的0所引起的。8.4.5 盲去卷积不以PSF为基础的图像复原方法统称为盲去卷积方法。其中一种常用的盲去卷积方法是以最大似然估计（MLE）为基础的，其最优解用规定的约束条件通过迭代并假设收敛来求解，得到的最大f(x,y)和h(x,y)就是还原的图像和PSF。（1）盲去卷积的实现工具箱通过函数deconvblind来实现：fr, PSFe=deconvblind(g, INITPSF)fr, PSFe=deconvblind(g, INITPSF, NUMIT, DAMPAR, WEIGHT) 其中INITPSF是点扩散函数的初始估计；PSFe是这个函数最终计算得到的估计值，PSF估计受其初始推测尺寸的影响巨大，而很少受其初始推测值的影响（一个值为1的数组是合理的初始推测）；fr是利用估计的PSF复原的图像，用来取得复原图像的算法是L-R迭代算法，默认迭代次数是10。(2)举例：使用deconvblind通过退化图像估计PSF/产生实际的退化PSFPSF=fspecial(gaussian, 7, 10);/产生一个大小为7*7、标准偏差为10的实际的退化PSF（高斯低通滤波模糊PSF）imshow(pixeldup(