朱正为8数字图像处理的应用.docx
《朱正为8数字图像处理的应用.docx》由会员分享，可在线阅读，更多相关《朱正为8数字图像处理的应用.docx（31页珍藏版）》请在冰豆网上搜索。
朱正为8数字图像处理的应用
第八章数字图像处理的应用
本章通过几个MATLAB数字图像处理实例，介绍数字图像处理技术的应用，具体内容包括：
☐MATLAB图像处理工具箱和基本的图像处理操作
☐图像直方图及其拟合
☐图像平滑：
中值滤波去噪
☐图像复原
8.1MATLAB图像处理工具箱和基本的图像处理操作
图像处理工具箱是一个函数的集合，它扩展了matlab数值计算环境的能力。
这个工具箱支持了大量图像处理操作，包括：
空间图像变换Spatialimagetransformations
形态操作Morphologicaloperations
邻域和块操作Neighborhoodandblockoperations
线性滤波和滤波器设计Linearfilteringandfilterdesign
格式变换Transforms
图像分析和增强Imageanalysisandenhancement
图像登记Imageregistration
清晰化处理Deblurring
兴趣区处理Regionofinterestoperations
工具箱里的函数都是M文件，可以通过typefunction_name来查看代码，也可以通过写自己的matlab函数来扩展工具箱。
要查看是否安装了该工具箱，可以使用ver命令来查看已安装的工具箱。
其他相关工具箱有：
DSPBlockset
ImageAcquisitionToolbox
MappingToolbox
SignalProcessingToolbox
WaveletToolbox
下面通过简单的例子介绍工具箱的使用方法，此例子介绍了一些基本的图像处理操作，包括读、写图像，演示图像的直方图均衡，得到图像信息等。
（1）MATLAB中的数字图像表示
在MATLAB中，一幅数字图像表示成矩阵（以二维数组的形式存储的）：
其中的每一个元素称为象素，符号f（p,q）表示位于第p行和q列的元素。
（2）读取一个图像
使用imread可以将图像读入MATLAB环境。
>>clear,closeall
>>f=imread（'D:
\myimages\dianluban.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）检查图像在工作区中的表现形式
>>whos
Name Size Bytes Class
I 291x240 69840 uint8array
Grandtotalis69840elementsusing69840bytes
图像数据是按数组的形式存储的。
（5）显示图像直方图和图像直方图均衡处理
显示图像的直方图
figure,imhist（I）
或figure,imhist（I,b）
//I是图像文件;b是灰度级的个数，省略默认为256
执行直方图均衡化
I2=histeq（I）;
figure,imshow（I2）
均衡后的直方图为
figure,imhist（I2）
（6）保存图像（写图像到硬盘）
>>imwrite（f,'filename','jpg'）
或imwrite（f,'C:
\pout2.png'）
（7）得到图像信息
>>imfinfo（'C:
\pout2.png'）//返回关于图像文件的信息
8.2图像直方图及其拟合
（1）概述
雷达目标和杂波幅度分布特性是雷达目标和杂波的主要统计特性之一，对雷达信号处理、检测、识别、仿真及对雷达系统设计和性能评估均有十分重要的意义。
长期以来，雷达工作者一直在研究和探讨这一问题，由于雷达杂波比较复杂，它包括地物杂波、海杂波、气象杂波、鸟群杂波（仙波）和箔条杂波等各种有源和无源干扰，并且在不同条件下，又是千变万化的，故分析起来比较困难。
一般都是采用统计的方法，对它们进行分析或对实测数据进行拟合来对雷达杂波的幅度分布进行建模。
到目前为止，雷达杂波幅度分布模型有高斯分布、瑞利分布、K分布、韦布尔分布、对数-正态分布、莱斯（Rice）分布模型等。
现代目标模型有莱斯模型、对数-正态分布模型和
（Chi方）模型等。
本节针对高分辨率SAR图像（的直方图），分析几种常用模型的建模能力（其核心是各分布概率密度函数参数的估计，其中常用的参数估计方法有矩估计法、极大似然估计法等）。
（2）需拟合的直方图
草地SAR图像、森林SAR图像和城市SAR图像。
（3）举例
采用对数正态分布对高分辨率SAR图像的直方图进行拟合。
（A）对数正态分布的定义
定义：
如果随机变量
的函数
服从正态分布
，
，则称
服从参数为
和
的对数正态分布，记作
～
。
对数正态分布的概率密度函数为：
（B）参数
与
的极大似然估计
根据极大似然估计法的原理，由上式，得似然函数为：
两边取对数，得：
似然方程组为：
故得
与
的极大似然估计分别是
其中
是
的无偏估计，
是
的有偏（渐进无偏）估计。
可以证明
的无偏估计为：
（C）MATLAB程序
（a）求对数正态分布参数的mu（
）和delta2（
）
clearall;
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;
%求mu
form=1:
p
forn=1:
q
ifx（m,n）==0
y（m,n）=0;
s=s+1;
else
y（m,n）=log（x（m,n））;
end
end
end
mu=sum（sum（y））/N;
%求sigma2
fori=1:
p
forj=1:
q
z（i,j）=（y（i,j）-mu）.^2;
end
end
sigma2=sum（sum（z））/N;
mu
sigma2
mu=4.8947,4.6667,3.8266
sigma2=0.0095,0.0341,0.0677
（b）绘制图像灰度直方图
clearall;
a=imread（'D:
\image002.jpg'）;
[p,q]=size（a）;
x=1:
1:
256;%设定X轴的范围
total=zeros（1,256）;%用于存放0～255各级灰度出现的次数。
zeros（m,n）将产生一个mxn的以零为项目之矩阵；而zeros（n）会产生一个nxn的以零为内容之方矩阵
b=double（a）;%将创建原图像的副本并将其中的数据元素转换为double型，uint8型的数据无法进行计算
sum=0;
fori=1:
p*q%统计0～255各级灰度在图像中出现次数total（b（i）+1）=total（b（i）+1）+1;%由于total的下标从1开始到256而b（i）的范围是0～255，所以需要将b（i）的值加1才可以使b（i）和total对应
end;
forj=1:
256%统计0～255各级灰度的累计出现次数
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（[xminxmaxyminymax]）
（c）绘制对数正态分布概率密度函数
clearall;
%mu=4.8947;
%sigma2=0.0095;
mu=4.6667;
sigma2=0.0341;
%mu=3.8266;
%sigma2=0.0677;
sigma=sqrt（sigma2）;
y=lognpdf（1:
300,mu,sigma）;
holdon;
plot（y,'b--'）;
xlabel（'灰度幅值'）;
ylabel（'概率密度'）;
8.3图像平滑：
中值滤波去噪
（1）在图像中加入椒盐噪声
>>f=imread（'E:
\pro7.jpg'）;%Theimagepro7.jpgislocatedintherootdirectoryofE.
>>type='salt&pepper';
>>p=0.2;
>>RN=imnoise（f,type,p）;
>>c=find（RN==1）;%将所有数组RN的索引返回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,[33],'symmetric'）;
>>imshow（h31）
>>h32=medfilt2（h31,[33],'symmetric'）;
>>figure,imshow（h32）
>>h33=medfilt2（h32,[33],'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=medfilt2（faddedgn,[77],'symmetric'）;
>>figure,imshow（hg1）
>>hg2=medfilt2（hg1,[77],'symmetric'）;
>>figure,imshow（hg2）
>>hg3=medfilt2（hg2,[77],'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模型化
●典型方法：
通过产生PSF及测试各种复原算法的结果来试验
●采用盲去卷积来推断PSF
本节讨论使用函数imfilter和fspecial，以及噪声生成函数imnoise建模PSF的技术。
由场景和传感器产生的图像模糊可用低通滤波器来建模；运动模糊使用IPT（图像处理工具箱）函数fspecial（生成线性空间滤波器）建模。
（1）产生运动模糊的PSF:
PSF=fspecial（'motion',len,theta）
调用fspecial将返回PSF，它近似于摄像机有len个像素的线性移动的效果。
参数theta以度为单位，以顺时针方向相对正水平轴进行度量。
（2）创建一个已知PSF的退化图像，使用函数imfilter（线性空间滤波）：
>>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次。
自定义函数pixeldup（A,m,n）的程序
functionB=pixeldup（A,m,n）
%PIXELDUPduplicatespixelsofanimageinbothdirections
%Checkinputs
ifnargin<2
error（'Atleasttwoinputsarerequired.'）;
end
ifnargin==2
n=m;
end
%Generatravectorwithelements1:
size（A,1）.
u=1:
size（A,1）;
%Duplicateeachelementofthevectormtimes.
m=round（m）;%Protectagainstnonintergers.
u=u（ones（1,m）,:
）;
u=u（:
）;
%Nowrepeatfortheotherdirection.
v=1:
size（A,2）;
%Duplicateeachelementofthevectormtimes.
n=round（n）;%Protectagainstnonintergers.
v=v（ones（1,n）,:
）;
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*512显示
8.4.2维纳滤波复原
维纳滤波复原是一种线性图像复原方法。
我们首先定义平均噪声功率和平均图像功率及它们的比率：
有时用
来近似代替
，以便产生一个常数数组，产生参数维纳滤波器，可以对直接逆滤波产生重大改进。
（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=deconvwnr（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（Sf）））;
Fr3=deconvwnr（g,PSF,NCORR,ICORR）;
结果如图（d）所示，所得结果虽然仍有一些噪声存在，但已和原图像很接近，因为原图像和噪声都是已知的，所以可以正确地估算参数。
在实践中，当这些量之一（或更多）未知时，挑战便是在实验中智能地选择所用的函数，直到获得可接受的结果为止。
8.4.3约束的最小二乘方滤波复原
约束最小二乘方滤波复原也是一种线性图像复原方法。
（1）约束最小二乘方滤波复原在IPT中通过函数deconvreg来实现：
fr=deconvreg（g,PSF,NOISEPOWER,RANGE）//NOISEPOWER与
成正比，其较好的初始估计为
，最后的结果可能有很大的不同。
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和RANGE进行实验，得出如图（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-R算法由函数deconvlucy实现：
fr=deconvlucy（g,PSF,NUMIT,DAMPAR,WEIGHT）
其中参数NUMIT为迭代次数（默认为10次）。
参数DAMPAR指定结果图像与原图像g之间的偏离阈值。
默认值为0。
当像素值偏离原值的范围在DAMPAR之内时，就不用再迭代。
WEIGHT是一个与输入图像g同样大小的数组（默认值是与g同样大小的一个单位数组），它为每一个像素分配一个权重来反映其逼真的程度。
例如一个不良像素会被赋以权重值零，从而paic该像素来求解。
当用一个指定的PSF来模拟模糊时，WEIGHT可从计算像素中剔除那些来自图像边界的像素点，若PSF的大小为n*n，则在WEIGHT中用到的零边界宽度为ceil（n/2）。
（2）举例
f=checkerboard（8）;//产生测试图像（a）
imshow（pixeldup（f,8））;//放大显示
PSF=fspecial（'gaussian',7,10）;//产生一个大小为7*7且标准偏差为10的高斯（低通滤波模糊）PSF
SD=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）/2）;/产生数组WEIGHT，WEIGHT数组的大小是64*64，且有值为0的4像素宽的边界，其余的像素都是1
WEIGHT=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次后的复原图像。
在所有结果图像中，细黑边界都是由数组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
//产生实际的退化PSF
PSF=fspecial（'gaussian',7,10）;//产生一个大小为7*7、标准偏差为10的实际的退化PSF（高斯低通滤波模糊PSF）
imshow（pixeldup（
升级会员 