1、第八章 数字图像处理的应用 本章通过几个 MATLAB 数字图像处理实例，介绍数字图像处理技术的应用，具体内容包括： MATLAB 图像处理工具箱和基本的图像处理操作 图像直方图及其拟合 图像平滑：中值滤波去噪 图像复原8.1 MATLAB 图像处理工具箱和基本的图像处理操作图像处理工具箱是一个函数的集合，它扩展了 matlab 数值计算环境的能力。这个工具箱支持了大量图像处理操作，包括：空间图像变换 Spatial image transformations形态操作 Morphological operations 邻域和块操作 Neighborhood and block operatio
2、ns线性滤波和滤波器设计 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 allf=imre
4、ad(D:myimagesdianluban.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
5、）检查图像在工作区中的表现形式 whosName Size Bytes ClassI 291x240 69840 uint8 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） 保存图像（
6、写图像到硬盘） imwrite(f, filename, jpg)或 imwrite(f, C:pout2.png)（7）得到图像信息imfinfo(C:pout2.png) / 返回关于图像文件的信息8.2 图像直方图及其拟合（1）概述雷达目标和杂波幅度分布特性是雷达目标和杂波的主要统计特性之一，对雷达信号处理、检测、识别、仿真及对雷达系统设计和性能评估均有十分重要的意义。长期以来，雷达工作者一直在研究和探讨这一问题，由于雷达杂波比较复杂，它包括地物杂波、海杂波、气象杂波、鸟群杂波（仙波）和箔条杂波等各种有源和无源干扰，并且在不同条件下，又是千变万化的，故分析起来比较困难。一般都是采用统计的
7、方法，对它们进行分析或对实测数据进行拟合来对雷达杂波的幅度分布进行建模。到目前为止，雷达杂波幅度分布模型有高斯分布、瑞利分布、K 分布、韦布尔分布、对数-正态分布、莱斯（Rice）分布模型等。现代 目标模型 有莱斯模型、对数-正态分布模型和 （Chi 方）模型等。本节针对高分辨率 SAR 图2像（的直方图） ，分析几种常用模型的建模能力（其核心是各分布概率密度函数参数的估计，其中常用的参数估计方法有矩估计法、极大似然估计法等）。（2）需拟合的直方图草地 SAR 图像、森林 SAR 图像和城市 SAR 图像。（3）举例采用对数正态分布对高分辨率 SAR 图像的直方图进行拟合 。（A）对数正态分布
8、的定义定义：如果随机变量 的函数 Xln服从正态分布 ，2,N，则称 服从参数为 和 的对数正态分布，记作 02X2X。,LN对数正态分布的概率密度函数为：（B ） 参数 与 的极大似然估计2根据极大似然估计法的原理，由上式，得似然函数为：两边取对数，得：似然方程组为：故得 与 的极大似然估计分别是2其中 是 的无偏估计， 是 的有偏（渐进无偏）估计。可以2证明 的无偏估计为：2（C）MATLAB 程序（a）求对数正态分布参数的 mu（ ）和 delta2（ ）2clear all;i1=imread(D:image003.jpg);i2=rgb2gray(i1); %将真彩色图像转换成灰度图
9、像x=double(i2); %创建原图像的副本并将其元素转换为double型，uint8型数据无法进行计算p=size(x,1);q=size(x,2);N=p*q;s=0;%求 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.8266sigm
10、a2=0.0095, 0.0341, 0.0677（b）绘制图像灰度直方图clear all;a=imread(D:image002.jpg);p,q=size(a);x=1:1: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
11、(b(i)+1)=total(b(i)+1)+1; %由于total的下标从1开始到256而b(i)的范围是0255，所以需要将b(i)的值加1才可以使b(i)和total对应end;for 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;
12、%sigma2=0.0095;mu=4.6667;sigma2=0.0341;%mu=3.8266;%sigma2=0.0677;sigma=sqrt(sigma2);y=lognpdf(1:300,mu,sigma);hold on;plot(y, b-);xlabel(灰度幅值);ylabel(概率密度);0 50 100 150 200 250 300 350 400 45000.0050.010.0150.020.0250.030.0350.04值值值值值值0 50 100 150 200 250 300 350 400 45000.0020.0040.0060.0080.010.01
13、20.0140.0160.0180.02值值值值值值0 50 100 150 200 250 300 350 400 45000.0050.010.0150.020.0250.030.0350.040.045值值值值值值8.3 图像平滑：中值滤波去噪（1）在图像中加入椒盐噪声 f=imread(E:pro7.jpg); % The image pro7.jpg is located in the root directory of E. type = salt p = 0.2; RN = imnoise(f, type, p); c = find(RN = 1); %将所有数组 RN 的索引返
14、回 c 中，它指向非零元素 fs = f; fs(c) =255; %在图像中加入盐噪声 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,s
15、ymmetric); imshow(h31) h32 = medfilt2(h31,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); imsho
16、w(f), figure, imshow(faddedgn) hg1 = medfilt2(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
17、-134 ）8.4.1 退化函数的建模成像设备实验法用数学方法把 PSF 模型化典型方法：通过产生 PSF 及测试各种复原算法的结果来试验采用盲去卷积来推断 PSF本节讨论使用函数 imfilter 和 fspecial，以及噪声生成函数imnoise 建模 PSF 的技术。由场景和传感器产生的图像模糊可用低通滤波器来建模；运动模糊使用 IPT（图像处理工具箱）函数 fspecial（生成线性空间滤波器）建模。（1）产生运动模糊的 PSF:PSF=fspecial(motion, len, theta)调用 fspecial 将返回 PSF，它近似于摄像机有 len 个像素的线性移动的效果。参
18、数 theta 以度为单位，以顺时针方向相对正水平轴进行度量。（2）创建一个已知 PSF 的退化图像，使用函数 imfilter（线性空间滤波）： g=imfilter(f, PSF, circular);其中 circular 为边界选项，用来减少边界效应。（3）使用函数 checkerboard 产生测试图案C=checkerboard(NP, M, N)其中 NP 是每个小正方形一边的像素数，M 是行数，N 是列数。测试板左半部分的亮正方形是白色的，右半部分的亮正方形是灰色的。若省略了 N，则默认为 M，若都省略，则默认为 8。（4）由于有些复原算法对于大图像很慢，因此通常对小图像做实验
19、，而通过像素复制以放大显示图像（自定义函数）：B=pixeldup(A, m, n)该函数将 A 的每个像素在垂直方向上总共复制了 m 次，在垂直方向上总共复制了 n 次。自定义函数 pixeldup(A,m,n)的程序function B=pixeldup(A,m,n)%PIXELDUP duplicates pixels of an image in both directions%Check inputsif nargin2error(At least two inputs are required.);endif nargin=2n=m;end%Generatr a vector wi
20、th elements 1:size(A,1).u=1:size(A,1);%Duplicate each element of the vector m times.m=round(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),
21、:);v=v(:);B=A(u,v);（5）举例：模糊噪声图像的建模f=checkerboard(8);/产生测试图像 (a)PSF=fspecial(motion, 7, 45);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
22、*512 显示8.4.2 维纳滤波复原维纳滤波复原是一种线性图像复原方法。我们首先定义平均噪声功率和平均图像功率及它们的比率： 1,AuvSvMN,Afuvf vARf有时用 来近似代替 ，以便产生一个常数数组，R,/,fSuv产生参数维纳滤波器，可以对直接逆滤波产生重大改进。（1）在 IPT 中，维纳滤波使用函数 deconvwnr(g, PSF)来实现：fr=deconvwnr(g, PSF) /假设噪信比为零，这种形式实现的即是逆滤波器。fr=deconvwnr(g, PSF, NSPR) /假设噪信功率比已知，这用于实现参数维纳滤波器。fr=deconvwnr(g, PSF, NACO
23、RR, FACORR) /假设噪声和未退化图像的自相关函数 NACORR 和 FACORR是已知的。由 可知，通过计2,fSuvFvfxyf算图像（或噪声）功率谱的傅里叶逆变换可得到图像（或噪声）的自相关函数 。,fxyf（2）举例：使用函数 deconvwnr 复原前述（运动）模糊噪声图像fr1=deconvwnr(g, PSF) /直接逆滤波，结果如图(b) 所示，这个结果由噪声的效果所决定。 fr2=deconvwnr(g, PSF, R) /使用比率 R 复原图像，结果如图(c)所示。比率 R 是利用噪声图像和原图像得到的，过程如下：Sn=abs(fft2(noise).2;nA=su
24、m(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(Sf);Fr3=deconvwnr(g, PSF, NCORR, ICORR);结果如图(d)所示，所得结果虽然仍有一些噪声存在，但已和原图像很接近，因为原图像和噪声都是已知的，所以可以正确地估算参数。在实践中，当这些量之一（或更多）未知时，挑战便是在实验中智能地选择所用的函数，
25、直到获得可接受的结果为止。8.4.3 约束的最小二乘方滤波复原约束最小二乘方滤波复原也是一种线性图像复原方法。（1）约束最小二乘方滤波复原在 IPT 中通过函数 deconvreg 来实现：fr=deconvreg(g, PSF, NOISEPOWER, RANGE)/ NOISEPOWER 与成正比，其较好的初始估计为 ，最后的结果可能有很2n 2MNm大的不同。RANGE 为值的范围，在求 时，该算法受这个范围的限制，默认范围是10 -9,109，MATLAB 中的符号为 1e-10, 1e10。若将上述两个参数排除在参量之外，则函数 deconvreg 就会产生一个逆滤波方案。（2）举例
26、：使用函数 deconvreg 复原前述（运动）模糊噪声图像因为图像大小为 64*64，噪声的方差为 0.001，均值为 0，所以NOISEPOWER 的初始估计为 。640.14图（A）显示了使用下述命令得到的结果：fr=deconvreg(g, PSF, 4);其中 g 和 PSF 来自 8.4.1 节例子。可以看出 NOISEPOWER 的值并不是非常好，通过对参数 NOISEPOWER 和 RANGE 进行实验，得出如图（B）所示的结果，它是使用如下命令得到的：fr=deconvreg(g, PSF, 0.4, 1e-7, 1e7);NOISEPOWER 的值下调了一个数量级，而 RA
27、NGE 比默认值更加紧缩。维纳滤波结果比该结果要好得多，但其条件是噪声和图像谱已知，若没有这些信息，则用两种滤波器通过实验得到的结果常常是可比的。8.4.4 使用 Lucy-Richardson 算法的迭代非线性复原（1）L-R 算法实现作为复原的工具，非线性迭代技术常常会获得比线性方法更好的结果。在 IPT 工具箱中，选择的非线性方法是由 Richardson 和 Lucy 开发的技术（L-R 算法是从最大似然公式引出来的） ，L-R 算法由函数deconvlucy 实现：fr=deconvlucy(g, PSF, NUMIT, DAMPAR, WEIGHT)其中参数 NUMIT 为迭代次数
28、（默认为 10 次） 。参数 DAMPAR 指定结果图像与原图像 g 之间的偏离阈值。默认值为 0。当像素值偏离原值的范围在 DAMPAR 之内时，就不用再迭代。WEIGHT 是一个与输入图像 g 同样大小的数组（默认值是与 g 同样大小的一个单位数组），它为每一个像素分配一个权重来反映其逼真的程度。例如一个不良像素会被赋以权重值零，从而 paic 该像素来求解。当用一个指定的PSF 来模拟模糊时，WEIGHT 可从计算像素中剔除那些来自图像边界的像素点，若 PSF 的大小为 n*n，则在 WEIGHT 中用到的零边界宽度为 ceil(n/2)。（2）举例f=checkerboard(8);/
29、产生测试图像 (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; /将 DAMPAR 赋以 10 倍于 SD 的值；LIM=ceil(size(PSF, 1)/
30、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 次后的复原图像。在所有结果图像中，细黑
31、边界都是由数组 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
32、是这个函数最终计算得到的估计值，PSF 估计受其初始推测尺寸的影响巨大，而很少受其初始推测值的影响（一个值为 1 的数组是合理的初始推测） ；fr是利用估计的 PSF 复原的图像，用来取得复原图像的算法是 L-R 迭代算法，默认迭代次数是 10。(2)举例：使用 deconvblind 通过退化图像估计 PSF/产生实际的退化 PSFPSF=fspecial(gaussian, 7, 10);/产生一个大小为 7*7、标准偏差为 10的实际的退化 PSF（高斯低通滤波模糊 PSF）imshow(pixeldup(PSF, 73), ); /放大显示 PSF/用 PSF 产生高斯低通滤波模糊和均
33、值为 0，标准偏差为 0.01 的高斯噪声污染的退化图像f=checkerboard(8);SD=0.01;g=imnoise(imfilter(f, PSF), gaussian, 0, SD.2);/使用 deconvblind 仅从给出的退化图像 g 获得 PSF 的估计值INITPSF=ones(size(PSF);NUMIT=5;/可同时取 10、20 比较一下结果DAMPAR=10*SD; LIM=ceil(size(PSF, 1)/2);WEIGHT=zeros(size(g);WEIGHT(LIM+1:end-LIM, LIM+1:end-LIM)=1;fr, PSFe=dec
34、onvblind(g, INITPSF, NUMIT, DAMPAR, WEIGHT);figure, imshow(pixeldup(PSFe, 73), );(a)原始 PSF；(b)到(d)在函数 deconvblind 中分别使用 5 次、10 次、20 次迭代估计得到的 PSF。1维纳滤波复原程序I=checkerboard(8);noise=0.1*randn(size(I);PSF=fspecial(motion,21,11);Blurred=imfilter(I,PSF,circular);BlurredNoisy=im2uint8(Blurred+noise);NP=abs(
35、fftn(noise).2;NPOW=sum(NP(:)/prod(size(noise);NCORR=fftshift(real(ifftn(NP);IP=abs(fftn(I).2;IPOW=sum(IP(:)/prod(size(noise);ICORR=fftshift(real(ifftn(IP);ICORR1=ICORR(:,ceil(size(I,1)/2);NSR=NPOW/IPOW;subplot(221);imshow(BlurredNoisy,);title(模糊和噪声图像);subplot(222);imshow(deconvwnr(BlurredNoisy,PSF,N
36、SR),);title(deconbwnr(A,PSF,NSR);subplot(223);imshow(deconvwnr(BlurredNoisy,PSF,NCORR,ICORR),);title(deconbwnr(A,PSF,NCORR,ICORR);subplot(224);imshow(deconvwnr(BlurredNoisy,PSF,NPOW,ICORR1),);title(deconbwnr(A,PSF,NPOW,ICORR_1_D);值值值值值值值 deconbwnr(A,PSF,NSR)deconbwnr(A,PSF,NCORR,ICORR) deconbwnr(A,P
37、SF,NPOW,ICORR1D)规则化滤波复原程序I=checkerboard(8);PSF=fspecial(gaussian,7,10);V=.01;BlurredNoisy=imnoise(imfilter(I,PSF),gaussian,0,V);NOISEPOWER=V*prod(size(I);J LAGRA=deconvreg(BlurredNoisy,PSF,NOISEPOWER);subplot(221);imshow(BlurredNoisy);title(A=Blurred and Noisy);subplot(222);imshow(J);title(J LAGRA=d
38、econvreg(A,PSF,NP);subplot(223);imshow(deconvreg(BlurredNoisy,PSF,LAGRA/10);title(deconvreg(A,PSF,0.1*LAGRA);subplot(224);imshow(deconvreg(BlurredNoisy,PSF,LAGRA*10);title(deconvreg(A,PSF,10*LAGRA);A=Blurred and Noisy J LAGRA=deconvreg(A,PSF,NP)deconvreg(A,PSF,0.1*LAGRA) deconvreg(A,PSF,10*LAGRA3Luc
39、y-Richardson 复原程序举例%Lucy-Richardson 复原程序I=checkerboard(8);PSF=fspecial(gaussian,7,10);V=.0001;BlurredNoisy=imnoise(imfilter(I,PSF),gaussian,0,V);WT=zeros(size(I);WT(5:end-4,5:end-4)=1;J1=deconvlucy(BlurredNoisy,PSF);J2=deconvlucy(BlurredNoisy,PSF,20,sqrt(V);J3=deconvlucy(BlurredNoisy,PSF,20,sqrt(V),
40、WT);subplot(221);imshow(BlurredNoisy);title(A=Blurred and Noisy);subplot(222);imshow(J1);title(deconvlucy(A,PSF);subplot(223);imshow(J2);title(deconvlucy(A,PSF,NI,DP);subplot(224);imshow(J3);title(deconvlucy(A,PSF,NI,DP,WT);A=Blurred and Noisy deconvlucy(A,PSF)deconvlucy(A,PSF,NI,DP) deconvlucy(A,PS
41、F,NI,DP,WT)4盲去卷积复原%盲去卷积复原I=checkerboard(8);PSF=fspecial(gaussian,7,10);V=.0001;BlurredNoisy=imnoise(imfilter(I,PSF),gaussian,0,V);WT=zeros(size(I);WT(5:end-4,5:end-4)=1;INITPSF=ones(size(PSF);FUN=inline(PSF+P1,PSF,P1);J P=deconvblind(BlurredNoisy,INITPSF,20,10*sqrt(V),WT,FUN,0);subplot(221);imshow(BlurredNoisy);title(A=Blurred and Noisy);subplot(222);imshow(PSF,);title(True PSF);
