05 第5章 MATLAB在数字图像处理中的应用.pptx

第1页：第5章MATLAB在数字图像处理中的应用

MATLAB在数字图像处理中有着广泛的用途。MATLAB提供的图像处理工具箱使得数字图像处理变得相当简单。以下几个例子包括图像配准、图像消噪、直方图均衡、图像压缩等图像处理的常用内容。希望通过这些例子，带领大家进入图像处理的基础领域。

第2页：本章知识架构

【本章知识架构】

第3页：基于FFT的图像配准·实验目的

1．实验目的

(1)进一步理解FFT的理论基础、仿真算法实现。

(2)了解常用的图像配准算法，并掌握各自的优缺点。

(3)掌握基于傅里叶变换的图像配准方法的基本原理及其仿真实验算法。2．实验原理由于图像的旋转、平移、比例变换都能在傅里叶变换频域中反映出来，利用变换域方法还有可能获得一定程度的抵抗噪声的鲁棒性，同时傅里叶变换可以采用FFT的方法提高执行的速度。

第4页：相位相关技术是配准两幅图像的平移失配的典

相位相关技术是配准两幅图像的平移失配的典型方法。相位相关方法依据傅里叶变换的平移性质。平移不影响傅立叶变换的幅值，无论图像沿x轴还是y轴平移，或者两个方向都平移，对应的幅值谱和原始图像都是完全一样的。

第5页：

第6页：

第7页：基于FFT的图像配准（续）

以上公式中的*表示共轭复数。由式(5-10)可见，由于两幅图像在离散空间域中存在着位移差，反映到频域中仅表现为相位上的差别。而这个相位差与图像间的平移量有直接关系。等式左边为一个虚指数,其傅里叶反变换恰好是一个狄拉克函数(空域中经过平移)，对其进行处理之后，其最大峰值点恰好位于x d和y d处。

第8页：基于FFT的图像配准·仿真思路

3．仿真思路

在图5.1中，x、y为配准前原始图像在水平和垂直方向上的平移量；x d、y d为通过计算之后待配准图像与原始图像在水平、垂直方向上存在的位移差值。对于本文中的配准方法，语言描述如下。

(1)配准前的图像准备，首先找到一幅参考图像，并分别对参考图像做3个方向上的平移(水平、垂直、与水平成45°角)。

(2)配准中的过程，读入两幅图像，分别为参考图像和待配准图像，调用配准程序算法，求出配准所需要的参数。

• (3)根据步骤
• (2)中得到的配准参数，对待配准图像进行平移，使得待配准图像恢复为原来的参考图像，从而达到图像的配准。

(4)参照理论数据，和实验数据进行比较，分析该种算法的性能。

第9页：基于FFT的图像配准·程序代码

4．程序代码

1)处理垂直方向配准的程序im1=imread(‘lena.bmp’); %读入lena图像，lena为256*256的图像im2=imread('lena7c.jpg'); %读入经过垂直移位7个像素后的lena图像subplot(1,2,1),imshow(im1);title('参考图像');subplot(1,2,2),imshow(im2);title('待配准图像');[xtan,ytran,delta]=computedelta(im1,im2); %调用配准函数，计算lena7c与lena的偏移a=imread('lena7c.jpg');for j=ytran:ytran+1;for x=1:256;for y=1:256;if (y-j)<0;k=y-j+256;elsek=y-j+1;endJ(k,x)=a(y,x);endendendfigure;imshow(J,[]);title('配准后图像');end

第10页：

第11页：处理水平方向配准的程序程序代码

2)处理水平方向配准的程序

im1=imread('lena.bmp'); %读入256*256的lena图像 im2=imread('lena7s.jpg'); %读入经过水平移位7个像素后的lena图像 subplot(1,2,1),imshow(im1);title('参考图像'); subplot(1,2,2),imshow(im2);title('待配准图像'); [xtan,ytran,delta]=computedelta(im1,im2); a=imread('lena7s.jpg'); for i=xtan:xtan+1; for y=1:256; for x=1:256; if (x-i)<0; k=x-i+256; else k=x-i+1; end J(y,k)=a(y,x); end end end figure; imshow(J,[]);title('配准后图像'); end

第12页：处理对角线方向配准的程序程序代码

3)处理对角线方向配准的程序

%用来处理对角线方向上的配准

function test() im1=imread('lena.bmp'); im2=imread('lena77.bmp'); subplot(1,2,1),imshow(im1);title('参考图像'); subplot(1,2,2),imshow(im2);title('待配准图像'); [xtan,ytran,delta]=computedelta(im1,im2); a=imread('lena77.bmp'); for i=xtan:xtan+1; for y=1:256; for x=1:256; if (x-i)<=0; k=x-i+256; else k=x-i+1; end J(y,k)=a(y,x);

第13页：基于FFT的图像配准·程序代码

end end end for j=ytran:ytran+1; for x=1:256; 256 for y=1:256; if (y-j)<0; k=y-j+256; else k=y-j+1; end J1(k,x)=J(y,x); end end end figure; imshow(J1,[]);title('配准后图像'); end

第14页：配准过程程序程序代码

4)配准过程程序

function [xtran,ytran,delta]=computedelta(im1,im2)

%估算两幅图像之间的平移差，返回水平方向平移差到xtran，返回垂直方向平移差到ytran

if(size(im1,3)==3) im1=rgb2gray(im1); im2=rgb2gray(im2); end im1=double(im1); im2=double(im2); F1=fft2(im1); F2=fft2(im2); fz=F1.*conj(F2); fm=abs(F1.*conj(F2)); div=fz./fm; delta=ifft2(div); figure mesh(delta);

第15页：基于FFT的图像配准·程序代码

axis([-10 10-10 10 0 1]); %通过求傅里叶反变换之后得到了一个狄拉克函数 maxvalue=max(max(delta)); %找出脉冲的最大值位置 delta=fliplr(flipud(delta)); %进行矩阵变换，便于计算最大脉冲所在位置 figure mesh(delta); axis([0 256-10 10 0 1]); %通过矩阵变换后的冲击函数 [ytran,xtran]=find(delta==maxvalue);

fprintf('图像2相对于图像1的位移量(单位-像素):\n水平方向(x方向) : %5.2f\n竖

直方向(y方向) : %5.2f\n',xtran,ytran); %实验结果保留5位有效数字 end

第16页：基于FFT的图像配准·运行结果与分析

5．运行结果与分析

• (1)将参考图像(图5.2)经垂直方向平移了7个整数像素点后，形成了如图5.3所示的待配准图像。

第17页：基于FFT的图像配准（续）

(2)对参考图像做水平方向7个整像素的平移，在这种情况下，形成的图像如图5.7和图5.8所示。从图5.8中看可以明显地看出它与图5.7在水平方向上存在着平移差。由程序仿真后求得平移量之差为x d=7.00，y d=0。这两个数据和先前做图像平移时的数据一致。并且由仿真后得到的配准后图像如图5.9所示。

第18页：

第19页：基于FFT的图像配准（续）

(3)将参考图像向与水平方成45°角方向平移7个单位。此时，得到的两幅图像如图5.10和图5.11所示。

第20页：基于FFT的图像配准（续）

通过实验仿真得到配准数据x d=7.00，y d =7.00，此数据正好是参考图像做平移的数据。两者一致则说明这两幅图像配准成功，配准后的图像如图5.12所示。通过此实验可以得出：在没有任何噪声的情况下，对于整数级像素点的傅里叶变换的配准有较高的配准精度。

第21页：基于FFT的图像配准·运行结果图

图5.12配准后图像

第22页：基于FFT的图像配准·思考题

6．思考题

(1)写出程序，要求

由一幅正常图像得到水平方向，垂直方向，水平、垂直方向皆有偏移的图像。

(2)给待配准图像加入噪声，用该方法进行配准，研究效果如何？如果效果不好，该如何改进？

第23页：基于场景的图像非均匀噪声抑·实验目的

1．实验目的

(1)了解红外探测器非均匀性的产生原因。

(2)学习基于场景的校正方法。

(3)掌握基于场景的校正方法中的恒定统计法。

第24页：基于场景的图像非均匀噪声抑·实验原理

2．实验原理造成图像非均匀性的因素是多方面的，起主要作用的是：回应率的非均匀性，包括光谱回应的非均匀性；读出电路自身及读出电路与探测器耦合的非均匀性；暗电流的非均匀性。图5.13标出了各种非均匀性的来源和表现形式。

第25页：基于场景的图像非均匀噪声抑·运行结果图

图5.13红外成像系统中各部分引入的非均匀性示意图

第26页：基于场景的图像非均匀噪声抑·仿真思路

3．仿真思路本例的结果是在MATLAB环境下用恒定统计法仿真所得的。本例中的图像序列是通过对一幅lena图像循环移位所产生的序列，在程序中将图像移位200次，将每100帧图像作为一个图像组，再对图像序列加噪声，由于乘性增益对图像的影响较小，所以在仿真过程中将Ak设为常数0.95，并将加性偏移Bk设为随机噪声，在MATLAB中由rand函数得到，由于固定噪声的缓变特性，所以可以将Bk认为在每一个图像组中是不变的。

第27页：

第28页：基于场景的图像非均匀噪声抑·程序代码

4．程序代码1)主程序

a = imread('lena.bmp'); [m,n]=size(a); b=round(10*rand(m,n)); HH=zeros(256,256);

%第1图像组

SNR1=[]; PSNR1=[]; for i=0:99 for x=1:256 for y=1:256 if (y+i)>256 k=y+i-256; else k=y+i; end J(x,k)=a(x,y); end end J=double(J); c=0.95; d=J*c+b; SNR1=[SNR1 snr(J,d)]; %计算校正前的100帧图像的信噪比总值

第29页：基于场景的图像非均匀噪声抑·程序代码

%PSNR1=[PSNR1 PSNR(J,d)]; d=uint8(d); pause(0.00001); subplot(1,2,1); imshow(d,[]); %加噪图像 %J=uint8(J); %fwname=sprintf('jiazaolena%2d.bmp',i); %imwrite(J,fwname); %保存图像 d=double(d); HH=HH+d; end SNR1=mean(SNR1); %计算平均信噪比 % PSNR1=mean(PSNR1); HH=HH/100; %计算输出信号平均值 ux=(0+255)/2; GG=ux*ones(256); B1=HH-c*GG; %估计噪声B % AB=B1-b; SNR2=[]; PSNR2=[]; for i=0:99 for x=1:256

第30页：基于场景的图像非均匀噪声抑·程序代码

for y=1:256 if (y+i)>256 k=y+i-256; else k=y+i; end J(x,k)=a(x,y); end end J=double(J); c=0.95; D=J*c+b-B1; % d=double(D); SNR2=[SNR2 snr(J,D)]; %计算校正后100帧图像的信噪比 %PSNR2=[PSNR2 PSNR(J,D)]; subplot(1,2,2); imshow(D,[]);title('校正后图像'); % D=uint8(D); % fwname=sprintf('quzaolena%2d.bmp',i); % imwrite(D,fwname); pause(0.0001); end

第31页：基于场景的图像非均匀噪声抑·程序代码

% AB=b-B1; %D=double(D); %R=RMSE(J,D); SNR2=mean(SNR2); %计算平均信噪比 %PSNR2=mean(PSNR2);2)仿真模型的均方根误差函数 function R=RMSE(T0,T1) tmp=0; for i=1:256 for j=1:256 tmp=tmp+(T1(i,j)-T0(i,j))*(T1(i,j)-T0(i,j)); end end tmp=sqrt(tmp); R=tmp/(256*256);

第32页：仿真模型的计算信号噪声比函数程序代码

3)仿真模型的计算信号噪声比函数

function snr=SNR(I,In) %I :original signal % In:noisy signal(ie. original signal +noise signal) % snr=10*log10(sigma2(I2)/sigma2(I2-I1)) [row,col]=size(I); snr=0; Ps=sum(sum(I.^2))/(256*256);%信号功率 Pn=sum(sum(((I-In).^2)))/(256*256);%噪声功率 snr=10*log10(Ps/Pn);

第33页：基于场景的图像非均匀噪声抑·运行结果与分

5．运行结果与分析第一组实验是对lena图像进行去噪，其结果如图5.15、图5.16和图5.17所示。其3幅图像结果是都是第50幅图像，加噪图像的信噪比是6.3404dB，去噪后图像的信噪比是10.8430dB，均方根误差RMSE位0.1236。从结果可以看出，恒定统计法可以对噪声图像进行去噪，且有一定的效果。

第34页：

第35页：第二组实验是对一幅均匀图像进行恒定统计法

第二组实验是对一幅均匀图像进行恒定统计法的去噪，其结果如图5.18、图5.19和图5.20所示。其结果都是第50幅图像，加噪图像的信噪比是8.0172dB，去噪后图像的信噪比是21.1590dB，均方根误差位0.0445。通过以上两组实验的对比可以看出，恒定统计法对均匀图像的去噪效果较好。

第36页：

第37页：基于场景的图像非均匀噪声抑·思考题

6．思考题

• (1)请列举该仿真的使用场景。
• (2)变换噪声的统计特性，重做此实验。

第38页：基于奇异值分解理论的图像去·实验目的

1．实验目的

• (1)学习奇异值分解的基本理论。
• (2)掌握基本的SVD算法，能够在MATLAB的环境下实现传统的基于多道信号的奇异值分解的去噪方法。
• (3)能使用MATLAB等语言实现基于Hankel矩阵的本征图像去噪方法。

第39页：

第40页：

第41页：

第42页：采用Hankel矩阵时

采用Hankel矩阵时，利用SVD可以将原始信号分解为一系列分量信号的一种简单线性叠加，这种简单线性叠加关系的优点是：一个分量从原始信号中被分离的过程就是从原信号中被简单地减去，这种减法运算将使得分离出来的各分量信号保持它们在原信号中的相位不变，即具有零相位偏移特性。式(5-23)也是采用Hankel矩阵时SVD的信号重构公式，它的意义还在于，可以选取感兴趣的几个分量进行简单的叠加，从而实现对信号特征信息的提取。

第43页：基于奇异值分解理论的图像去·仿真思路

3．仿真思路奇异值分解去噪基于矩阵奇异值分解的理论，将一个多道地震记录看成一个图，并将其分解为由一系列正交向量组成的子图，通过方差分析法或相干法，对子图进行选择并重构，达到去噪的目的。奇异值分解去噪利用了信号与噪声的能量可分性。在处理地震资料时，待处理区域的各道信号相关性较高，经过SVD(奇异值分解)后，信号能量越集中，信号与噪声越易分离；反之，信号与噪声不易分离。常规的SVD滤波只能提高水平同相轴的信噪比，而不能有效提高非水平或弯曲同相轴的信噪比，而如果采用多项式拟合的方法，检测出信号同相轴多项式，然后按照多项式确定的同相轴形状重新选取数据矩阵，再进行SVD滤波。这是对常规方法的改进。它的去噪手段灵活、保真性好、分辨率高，对一些突然的脉冲干扰、侧反射以及其他的反向干扰有明显效果。

第44页：基于奇异值分解理论的图像去·程序代码

4．程序代码

1) SVD去噪法的程序代码

clc clear all close all N=300; L=0.1; detat=L/N; t=(0:N-1)*detat; signal=Ricker(40,t-0.05); seismic=zeros(2*N,N/4);for k=1:floor(N/4) seismic(1:N,k)=signal;

第45页：基于奇异值分解理论的图像去·程序代码

end for k=1:floor(N/4) seismic(4*k+1:4*k+N,k)=seismic(4*k+1:4*k+N,k)+signal'; end figure(1) wigb(seismic) xlabel('trace') ylabel('time/ms') title('理想信号') [m n]=size(seismic); complex=seismic+0.155*randn(m,n);//加了0.155倍的随机噪声 figure(2) wigb(complex) xlabel('trace') ylabel('time/ms') title('加噪信号') [U,S,V]=svd(complex);figure(3)

S1=diag(S);

第46页：基于奇异值分解理论的图像去·程序代码

stem(S1); xlabel('V1i') ylabel('value') title('奇异值的大小') M=input('分界点：'); for i=1:75 if i>M

S(i,i)=0;end seismic2=U*S*V'; figure(4); wigb(seismic2) xlabel('trace') ylabel('time/ms') title('去噪信号')

第47页：Hankel矩阵去噪法的代码程序代码

2) Hankel矩阵去噪法的代码

clc clear all close all N=300; L=0.1; detat=L/N; t=[0:N-1]*detat; signal=Ricker(40,t-0.05); seismic=zeros(2*N,N/4); for k=1:floor(N/4) seismic(1:N,k)=signal; end for k=1:floor(N/4) seismic(4*k+1:4*k+N,k)=seismic(4*k+1:4*k+N,k)+signal'; end

第48页：基于奇异值分解理论的图像去·程序代码

figure(1) wigb(seismic) title('理想信号') [m n]=size(seismic); complex=seismic+0.155*randn(m,n); figure(2) wigb(complex) title('加噪信号') for i=1:75 y=complex(:,i); %取出每一列信号单独做处理 yout=cadzow1(y,2); %将单列信号排成一个Hankel矩阵，在cadzow函数中处理

seismic2(:,i)=yout; %恢复出进行过SVD去噪的单列信号，排在原来的位置上

end figure(3); wigb(seismic2) title('去噪信号')

第49页：基于奇异值分解理论的图像去·运行结果与分

5．运行结果与分析加了随机噪声后的信号图如图5.21所示。图5.21加了随机噪声后的信号图

第50页：用SVD法去噪之后的地震信号图

用SVD法去噪之后的地震信号图如图5.22所示。图5.22用SVD法去噪之后的地震信号图

第51页：Hankel矩阵法去噪后的信号图

Hankel矩阵法去噪后的信号图如图5.23所示。图5.23 Hankel矩阵法去噪后的信号图

第52页：基于奇异值分解理论的图像去·思考题

6．思考题

• (1)奇异值分解去噪法的原理是什么？
• (2)给lena图加随机噪声，用奇异值分解法去噪，并计算去噪前后的信噪比。

第53页：基于小波理论的图像消噪处理·实验目的

1．实验目的

(1)掌握小波变换的基本理论。

(2)研究信号去噪技术。

(3)进行计算机的仿真试验，提供试验程序及结果。

第54页：基于小波理论的图像消噪处理·实验原理

2．实验原理1)小波变换概述小波变换是一种信号的时间——尺度分析方法，它具有多分辨率分析的特点，在时频两域都具有表征信号局部特征的能力，是一种窗口大小固定不变但其形状可变，时间窗和频率窗都可以改变的时频局部化分析方法。在低频部分具有较高的频率分辨率，在高频部分具有较高的时间分辨率和较低的频率分辨率，很适合于探测正常信号中夹带的瞬态反常现象并展示其成分，所以被誉为分析信号的显微镜，利用连续小波变换进行动态系统故障检测与诊断具有良好的效果。

第55页：

第56页：基于小波理论的图像消噪处理·仿真思路

3．仿真思路本例研究小波阈值去噪，使用MATLAB进行仿真，使用的图像是256*256灰度图Lena，加入的噪声为高斯白噪声，用小波阈值法进行去噪处理。峰值信噪比(PSNR)是图像去噪效果的一个重要的参数指标，对于尺寸为N*N，量化级为0: 255的图像的PSNR定义为：

第57页：基于小波理论的图像消噪处理（续）

其中f′(x,y)为重建后的图像，峰值信噪比越高，则说明图像的去噪效果越好。此峰值信噪比公式适用于灰度图像。

第58页：基于小波理论的图像消噪处理·程序代码

4．程序代码

%读入原始图像并显示

i=imread('C:\\lena.jpg'); %imread函数的参数为图像所在绝对路径 figure(1)

%画出原始图像

subplot(2,2,1); imshow(i); title('原始图像'); axis square;

%生成含噪图像并显示

j=imnoise(i,'gaussian',0,0.01);

%画出加噪图像

subplot(2,2,2); imshow(j); title('含噪图像'); axis square;

%用sym4小波函数对j进行2层分解

[c,l]=wavedec2(j,2,'sym4');

第59页：基于小波理论的图像消噪处理·程序代码

%实现低通滤波消噪

a1=uint8(wrcoef2('a',c,l,'sym4',2));

%用coif2小波函数对低频系数a1进行2层分解

[gc,gl]=wavedec2(a1,2,'coif2'); n=[1,2];%设置尺度向量 p=[10.28,24.08];%设置阈值向量

%对3个高频系数进行阈值处理

nc=wthcoef2('h',gc,gl,n,p,'s'); nc=wthcoef2('v',gc,gl,n,p,'s'); nc=wthcoef2('d',gc,gl,n,p,'s'); mc=wthcoef2('h',gc,gl,n,p,'s'); mc=wthcoef2('v',gc,gl,n,p,'s'); mc=wthcoef2('d',gc,gl,n,p,'s');

%对阈值处理的小波系数进行重建

x2=waverec2(mc,gl,'coif2');

%显示去噪图像

subplot(2,2,3); imshow(uint8(x2)); title('消噪图像'); axis square;

第60页：基于小波理论的图像消噪处理·程序代码

%计算psnr q=0.0; for x=1:256; for j=1:256; diff=x2(x,j)-double(i(x,j)); q=q+diff*diff; end end psnr_y = q / (256*256); psnr_y = (255 * 255)/psnr_y; psnr=10*log10(psnr_y)

第61页：基于小波理论的图像消噪处理·运行结果与分

5．运行结果与分析小波阈值法可以很好地去除高斯白噪声，降噪图像具有高信噪比，其结果如图5.25所示。上例中得到消噪后图像的PSNR见表5-1。

第62页：基于小波理论的图像消噪处理·思考题

6．思考题

• (1)小波基的选择对去噪效果有何影响？
• (2)小波分解的级数对去噪效果有何影响？

第63页：彩色图像灰度化并进行直方图·实验目的

1．实验目的

(1)了解直方图和直方图均衡的定义。

(2)学会对图像进行灰度化、均衡化。

(3)学会利用灰度变换法提高对比度。

第64页：

第65页：彩色图像灰度化并进行直方图·仿真思路

3．仿真思路

(1)通过对一个256*256的lena彩色图像进行灰度化处理。

(2)绘制灰度图像的直方图(假设灰度级k=256)。

(3)利用灰度图的直方图绘制均衡直方图和灰度变化曲线。

(4)画出均衡化图像。

第66页：彩色图像灰度化并进行直方图·程序代码

4．程序代码

%输入文件

lena.jpg待处理图像

%输出图形窗口说明

% figure NO 1待处理彩色图像 % figure NO 2灰度化后图像 % figure NO 3直方图 % figure NO 4均衡化后直方图 % figure NO 5灰度变化曲线 % figure NO 6均衡化后图像 clear all

第67页：图像的预处理程序代码

1)图像的预处理

PS=imread(‘lena.jpg’); %读入JPG彩色图像文件 figure,imshow(PS) %显示出来figure NO 1 title('输入的彩色JPG图像') PW=rgb2gray(PS); %灰度化后的数据存入数组

2)绘制直方图

[m,n]=size(PW); %测量图像尺寸参数 GP=zeros(1,256); %预创建存放灰度出现概率的向量 for k=0:255 GP(k+1)=length(find(PW==k))/(m*n); %计算每级灰度出现的概率，将其存入GP中相应位置 end figure, subplot(2,1,1) bar(0:255,GP,'g') %绘制直方图figure NO 3 title('原图像直方图') xlabel('灰度值') ylabel('出现概率')

第68页：直方图均衡化程序代码

3)直方图均衡化

S1=zeros(1,256); for i=1:256 for j=1:i S1(i)=GP(j)+S1(i); %计算Sk end end S2=round(S1*256); %将Sk归到相近级的灰度 for i=1:256 GPeq(i)=sum(GP(find(S2==i))); %计算现有每个灰度级出现的概率 end subplot(2,1,2), bar(0:255,GPeq,'b') %显示均衡化后的直方图figure NO 4

title('均衡化后的直方图')

xlabel('灰度值') ylabel('出现概率') figure,plot(0:255,S2,'r') %显示灰度变化曲线figure NO 5 legend('灰度变化曲线') xlabel('原图像灰度级') ylabel('均衡化后灰度级')

第69页：图像均衡化程序代码

4)图像均衡化

PA=PW; for i=0:255 PA(find(PW==i))=S2(i+1); %将各个像素归一化后的灰度值赋给这个像素 end figure, subplot(2,1,1),imshow(PW) %显示灰度化后的图像，也是均衡化前的样品figure NO 2 title('灰度化后的图像') subplot(2,1,2),imshow(PA) %显示均衡化后的图像figure NO 6 title('均衡化后图像')

第70页：彩色图像灰度化并进行直方图·运行结果与分

5．运行结果与分析程序输出的结果如图5.26、图5.27、图5.28、图5.29、图5.30和图5.31所示。图5.26 lena原图

第71页：彩色图像灰度化并进行直方图·运行结果图

图5.27原图像的灰度级

第72页：彩色图像灰度化并进行直方图·运行结果图

图5.28原图像直方图

第73页：彩色图像灰度化并进行直方图·运行结果图

图5.29均衡化后的直方图

第74页：原图像灰度化后

原图像灰度化后，图像的直方图和均衡化的直方图对比，原灰度图像的直方图变换成均匀分布的形式。图像灰度化后，原图像与均衡化的图像相比，增加了图像灰度值的动态范围，从而提高了图像的整体对比度，使得均衡后的图像更加清晰。图5.30灰度化后的图像图5.31均衡化后的图像

第75页：彩色图像灰度化并进行直方图·思考题

6．思考题图像均衡的方法除了直方图均衡外，还有什么别的方法？试用MATLAB实现。

第76页：彩色电视图像生成原理·实验目的

1．实验目的

(1)了解电视扫描原理。

(2)了解电视机中三基色基本原理。

(3)学会用MATLAB软件进行电视扫描原理的仿真和彩色电视图像的表示。

第77页：彩色电视图像生成原理·实验原理

2．实验原理1)电视隔行扫描原理在电子显像管上，水平扫描电压和垂直扫描电压控制显像管电子枪发出的电子束轰击到显示屏的位置，而电子枪电压则控制电子束轰击荧光粉的速度，从而控制荧光粉发光亮度。电子束从左到右，从上到下周期地逐一快速扫过显示屏，利用荧光粉发光的暂时记忆效应形成图像。在模拟电视标准中，为了节约传输频带，采用了隔行扫描技术，即把一幅图像分为两次扫描显示，第一次扫描显示奇数行，第二次扫描显示偶数行。

第78页：三基色原理

2)彩色电视三基色原理白光通过棱镜后被分解成多种颜色逐渐过渡的色谱，颜色依次为红、橙、黄、绿、青、蓝、紫，这就是可见光谱。其中人眼对红、绿、蓝最为敏感，人的眼睛就像一个三色接收器的体系，大多数的颜色可以通过红、绿、蓝三色按照不同的比例合成产生。同样绝大多数单色光也可以分解成红绿蓝3种色光。这是色度学的最基本原理，即三基色原理图5.32。3种基色是相互独立的，任何一种基色都不能由其他两种颜色合成。红绿蓝是三基色，这3种颜色合成的颜色范围最为广泛。红绿蓝三基色按照不同的比例相加合成混色称为相加混色。

例：相加合成混色称为相加混色。

红色+绿色=黄色

绿色+蓝色=青色

红色+蓝色=紫色

红色+绿色+蓝色=白色图5.32三基色原理

第79页：彩色电视图像生成原理·仿真思路

3．仿真思路在数学上，图像可以用矩阵来表示，矩阵中的元素表示图像的一个点，称为像素，图像的总像素等于矩阵中元素的总个数。对于灰度图像，矩阵元素的值等于对应元素的亮度电平值。通常规定的亮度取值范围为0～255之间，用来表示256级灰度，0表示最暗，255表示最亮，这样一个像素就可以用8bit来表达。MATLAB中，用函数imread来读取图像文件，imshow将读取的图像矩阵显示到作图窗口，而imwrite可将图像矩阵保存为指定格式的图像文件。为了节约内存，这3个函数中的图像矩阵元素采用8bit无符号整型数(uint8)表示。

第80页：彩色电视图像生成原理·程序代码

4．程序代码

1) 8灰度条图像：eightgraybar.m %eightgraybar.m Y=[1:-1/7:0].*255; %8级灰度 N=8; %8*N列数=8*8=64列 S=ones(N,1); Y=S*Y; Y=reshape(Y,1,N*8); %按列优先提取Y中1*64个元素，返回1*64结构的Y矩阵， M=40; %M行数为40行 S=ones(M,1);

Y=S*Y; %重复M次矩阵行，得到图像矩阵

I=uint8(Y); %转换为uint8(无符号整型数据)格式

%显示并保存为图像文件

imshow(I); imwrite(I,'graybar.tif','tif');

第81页：隔行扫描程序代码

2)隔行扫描：gehangsaomiao.m % gehangsaomiao.m I=imread('graybar.tif'); %读入图像 img=double(I')/255; %转换为双精度数以便计算 [M,N]=size(img); for K=1:2 %K=1扫描奇数行；K=2扫描偶数行 for column=K:2:N for row=1:M dot=img(row,column); %取当前扫描点的灰度值 greydot=[dot,dot,dot]; %R、G、B设置为相同值，即显示灰度 plot(row,N-column,'s','MarkerFaceColor',greydot, 'MarkerEdgeColor',greydot,'Markersize',6); hold on; %保持扫描点 axis([1,M,1,N]); %设置作图范围 set(gcf,'DoubleBuffer','on'); %双缓冲避免作图闪烁 drawnow; %立即作图 end end end

第82页：彩条图像程序代码

3)彩条图像：Colorbar.m %Colorbar.m:生成400*480彩条图像 clear; R=[1 1 0 0 1 1 0 0]; G=[1 1 1 1 0 0 0 0]; B=[1 0 1 0 1 0 1 0]; % N*8列 N=60; S=ones(N,1); R=S*R; R=reshape(R,1,N*8);%按列优先提取R中1*480个元素，返回1*480结构的R

矩阵

G=S*G; G=reshape(G,1,N*8);%按列优先提取G中1*480个元素,返回1*480结构

的G矩阵

B=S*B; B=reshape(B,1,N*8);%按列优先提取B中1*480个元素，返回1*480结构的B

矩阵

M=400; %行数 S=ones(M,1);

%得出对应区域的三基色矩阵

第83页：彩色电视图像生成原理·程序代码

R=S*R; G=S*G; B=S*B;

%组合为三维数组

I(:,:,1)=R; I(:,:,2)=G; I(:,:,3)=B;

%类型转换

将I转换为8bit无符号整型数(uint8)数据类型显示出来

I=uint8(I*255);

%显示并保存为图像文件

imshow(I); imwrite(I,'colorbar.tif','tif' );

第84页：彩色电视图像生成原理·运行结果与分析

5．运行结果与分析8级灰条图如图5.33所示。图5.33 8级灰条图

第85页：彩色电视图像生成原理·运行结果图

图5.34为电视隔行扫描过程仿真得到的图像。仿真结果中动态的扫描过程符合电视机隔行扫描理论分析的扫描过程。即先扫描奇数行，后扫描偶数行，从而完成一幅图像的扫描过程。图5.34隔行扫描仿真：奇数行扫描完毕，正在扫描偶数行

第86页：彩色电视图像生成原理·运行结果图

图5.35为利用彩色电视的三基色原理仿真电视图像中经典的ColorBar彩条图像，仿真结果与理论的三基色混色原理完全符合。即红色+绿色=黄色；绿色+蓝色=青色；红色+蓝色=紫色；红色+绿色+蓝色=白色。图5.35彩条图像的产生

第87页：彩色电视图像生成原理·运行结果图

图5.35生成400*480像素的彩条图像，其彩色图像的一个像素由红(R)、绿(G)、蓝(B)3种颜色分量组合而成。相应地，表示像素的矩阵元素就是一个三元素组成的向量(R,G,B)。6．思考题

• (1)彩色图像在MATLAB中是如何表示的？
• (2)读入lena的彩色图像，分别分离出其中的红色、绿色、蓝色部分，并分别显示出来。

第88页：基于小波变换的图像压缩编码·实验目的

1．实验目的

(1)了解小波变换在图像压缩中的应用原理。

(2)了解MATBLAB中的图像处理工具箱及小波变换工具箱的使用。

第89页：基于小波变换的图像压缩编码·实验原理

2．实验原理小波变换用于图像编码的基本思想是：把图像进行多分辨率分解，分解成不同空间、不同频率的子图像，然后对子图像进行系数编码。和变换编码相同，基于小波变换的编码主要分为两步：

• ①小波变换；
• ②变换后系数处理。图像经过小波变换后生成的小波图像数据总量与原图像数据总量相等，即小波变换本身不具有压缩功能。之所以用于图像压缩，是因为小波变换后系数主要集中在低频部分，而水平、垂直和对角线部分能量较少。小波变换后系数处理方法一般为量化、重排列、熵编码。量化是根据肉眼视觉系统特性，减少肉眼不敏感的高频成分以达到压缩数据的目的，这一步骤对图像变换有损。重排列是根据频带分布相似性或重要性级别等性质将数据重新排列，是一种减少数据量的操作，为下一步的熵编码提供了很好的数据特性。熵编码是为了进一步减少数据量而采取的一种无损压缩操作。

第90页：基于小波变换的图像压缩编码·仿真思路

3．仿真思路将小波分析引入图像压缩的范畴也是一个重要的手段，并且它有着自己的特点。它的特点在于压缩比高、压缩速度快，压缩后能保持信号与图像的特征基本不变，且在压缩过程中可以抗干扰等。一个图像做小波分解后，可得到一系列不同分辨率的子图像，各个子图像对应的频率是不相同的。高分辨率(高频)子图像上大部分点的数值都接近于0，越是高就越是明显。而对于一个图像来说，表现一个图像的最主要的部分是低频部分，所以最简单的压缩方法是利用小波分解去掉图像的高频部分而只保留低频部分。

第91页：基于小波变换的图像压缩编码·程序代码

4．程序代码

load wbarb;

%装入图像

subplot(221); image(X);

%显示图像

colormap(map) title('原始图像'); axis square disp('压缩前图像X的大小'); whos('X') [c,s]=wavedec2(X,2,'bior3.7');

%对图像用小波进行层分解

ca1=appcoef2(c,s,'bior3.7',1);

%提取小波分解结构中的低频系数和高频系数

ch1=detcoef2('h',c,s,1);

%水平方向

cv1=detcoef2('v',c,s,1);

%垂直方向

cd1=detcoef2('d',c,s,1);

%斜线方向

第92页：基于小波变换的图像压缩编码·程序代码

a1=wrcoef2('a',c,s,'bior3.7',1); h1=wrcoef2('h',c,s,'bior3.7',1); v1=wrcoef2('v',c,s,'bior3.7',1); d1=wrcoef2('d',c,s,'bior3.7',1);

%各频率成分重构

c1=[a1,h1;v1,d1]; subplot(222); image(c1);

%显示分频信息

axis square;

title('分解后低频和高频信息');

%进行图像压缩

%保留小波分解第一层低频信息

%首先对第一层信息进行量化编码

ca1=appcoef2(c,s,'bior3.7',1); ca1=wcodemat(ca1,440,'mat',0); ca1=0.5*ca1; subplot(223); image(ca1);

%改变图像高度并显示

colormap(map);

第93页：基于小波变换的图像压缩编码·程序代码

axis square; title('第一次压缩图像');

disp('第一次压缩图像的大小为：');

whos('ca1') ca2=appcoef2(c,s,'bior3.7',2);

%保留小波分解第二层低频信息进行压缩

ca2=wcodemat(ca2,440,'mat',0);

%首先对第二层信息进行量化编码

ca2=0.25*ca2;

%改变图像高度并显示

subplot(224); image(ca2); colormap(map); axis square; title('第二次压缩图像');

disp('第二次压缩图像的大小为：');

whos('ca2')

第94页：

第95页：基于小波变换的图像压缩编码（续）

在这里可以看出，第一次压缩是提取原始图像中小波分解第一层的低频信息，此时压缩效果比较好，压缩比较小(约为1/3大小)。第二次压缩是提取第一层分解低频部分的低频部分(即第二层的低频部分)，其压缩比较大(约为1/2大小)，压缩效果在视觉效果上也基本过得去。上面的保留原始图像中低频信息的压缩办法只是一种最简单的压缩办法。它不需经过其他处理即可获得较好的压缩效果。当然，对于上面的实验还可以只提取小波分解的第三层和第四层的低频信息。从理论上说，可以获得任意压缩比的压缩图像。只不过在对压缩比和图像质量都有较高要求时，它就不如其他的编码方法了。用小波变换对图像进行压缩的结果如图5.36所示。

第96页：基于小波变换的图像压缩编码·运行结果图

图5.36用小波变换对图像进行压缩

第97页：基于小波变换的图像压缩编码·思考题

6．思考题

(1)利用小波变换压缩图像是有损压缩还是无损压缩？

(2)利用小波变换压缩图像主要利用了图像的什么特点？

第98页：美丽的分形图案·实验目的

1．实验目的

(1)了解分形问题的基本概念及常识。

(2)知道分形问题数学上是如何产生的。

(3)对于给出算法的分形问题，能够用计算机实现分形图形的绘制。2．实验原理在20世纪60年代，B.Mandelbrot对一个具有复杂几何性质但局部看起来仍然一样的几何对象提出了分形的概念。今天，在非线性动力学动力系统等许多领域都会看到分形的例子。随着电子计算机的发展，许多美丽的分形图案被构造出来，这个实验主要介绍两个分形图案。下面给出分形图案中的数据计算方法，范例1中的蕨类植物生长模型采用仿射变换的循环迭代方法。设x是一个n维向量，A是(n，n)维矩阵，b是与A同维的向量，则x→Ax+b在数学上称为仿射变换，取不同的A、b就会得到不同的变换结果。范例2中的图像同样是采用循环的仿射变换绘制的，只是参数选择不同而已。

第99页：美丽的分形图案·程序代码

3．程序代码

1)蕨类植物模拟

x = [.5; .5]; %初值 h = plot(x(1),x(2),'.'); %绘制初值点

%设置用于后面随机数的判别向量

p = [ .85 .92 .99 1.00]; b1 = [0; 1.6]; b2 = [0; 1.6]; b3 = [0; .44];

%仿射变换矩阵

A1 = [ .85 .04;-.04 .85]; A2 = [ .20-.26; .23 .22]; A3 = [-.15 .28; .26 .24]; A4 = [ 0 0 ; 0 .16]; for i=1:20000 r = rand; %产生随机数 if r < p(1) x = A1*x + b1; elseif r < p(2) x = A2*x + b2;

第100页：美丽的分形图案·程序代码

elseif r < p(3) x = A3*x + b3; else x = A4*x; end plot(x(1),x(2),'g'),hold on %采用绿色绘制 end

2)编写Sierpinski三角形分形图案函数

%Sierpinski三角形分形图案

%初值

x = [0; 0]; h = plot(x(1),x(2),'.'); %绘制初值 A=[1/2 0;0 1/2]; %仿射变换矩阵 b1=[0 0]'; b2=[1/2 0]'; b3=[1/4 sqrt(3)/4]'; for i=1:10000 r = rand; %产生随机数 if r >0.66666666667 x = A*x + b1; elseif r >0.333333333333

第101页：美丽的分形图案·程序代码

x=A*x+b2; else x = A*x+b3; end plot(x(1),x(2),'-'),hold on %保持图形 end axis off %取消坐标轴

4．运行结果与分析

蕨类植物的生长模拟图形如图5.37所示。图5.37蕨类植物的生长模拟图形

第102页：Sierpinski三角形程序代码

Sierpinski三角形如图5.38所示。图5.38 Sierpinski三角形

第103页：美丽的分形图案（续）

这个实验重点在于介绍如何绘制分形图案的方法，对于分形中的复杂数学问题没有进行深入的讨论，这些内容当然也不是仅仅一个实验可以解决的。但通过这两个分形图像的绘制，给出了典型的分形图案绘制方法，希望可以对以后的分形图案绘制做些铺垫。5．思考题

• (1)构成一副分形图像最基本的内容是什么？
• (2)变换不同的旋转仿射系数，重构不同的分形图像。

第104页：分形图像压缩编码·实验目的

1．实验目的

(1)了解分形图像压缩编码的原理及实现方法。

(2)了解简库恩算法。

(3)会利用MATLAB进行分形图像压缩编码。2．实验原理在现实生活中，图像是非严重自相似的，即严格分形(具有整体与局部自相似性)图像在现实生活中只占极小部分。因此，Barnsly提出将原图像预分解为若干个分形子块，使每一个子块具有一定的分形结构，即子块的整体与局部间存在某种自仿射特征的思路。但是他没有解决：用什么方法分割，在图像分割时，计算机怎样自动分割，分形库如何建等问题。因此，用他提出的方法是很难实现自动压缩的。在此基础上，简恩提出了基于图像分块方式来实现以局部的仿射变换代替全局的仿射变换的全自动分形编码方法，这是目前分形图像压缩编码中被主要使用的方法。

第105页：分形图像压缩编码·程序代码

4．程序代码1)分形图像压缩编码

lena=imread('lena1.tif'); I=double(lena);

N=256; %整个图像每行的像素

M=256;

J=4; %子块每行的像素个数

D=2*J; %8; %父块每行的像素 r=ones(N/J,M/J); %记录对应于每个R(i，j)的r值 h=zeros(N/J,M/J); %记录对应于每个R(i，j)的h值 nn=ones(N/J,M/J); %记录对应于每个R(i，j)的仿射变换号n DD1=ones(N/J,M/J); %记录对应于每个R(i，j)的父块D左上角的行号 DD2=ones(N/J,M/J); %记录对应于每个R(i，j)的父块D左上角的列号 R=zeros(J,J); B=zeros(D,D); B1=zeros(J,J); B2=zeros(J,J); for i=1:J:(N-J+1)

for j=1:J:(M-J+1) %第一层循环的目的是遍历所有的子块

第106页：分形图像压缩编码·程序代码

e=100; R=out(I,J,i,j); %用out函数从原图像I中取出一个子块R，i,j为子块R的左上

角坐标

x=floor((i-1)/J)+1; y=floor((j-1)/J)+1; for k=1:(N-D+1) %test %第二层循环是针对每一子块R(i,j)为寻找与之相仿的父块 for l=1:(M-D+1) % test%k ,l为父块的左上角坐标 B=out(I,D,k,l); %取块操作

B1=average(B,D); %把大小为D*D的父块通过取平均像素变为大小为J*J的子

块B1 for n=1:8 %对B1进行8种仿射变换之一得到B2 B2=Ln(B1,J,n); %对B1进行8种仿射变换，得到B2 [r0,h0,dp0]=rhdp(R,B2,J); if dp0<e e=dp0; r(x,y)=r0; h(x,y)=h0; nn(x,y)=n;

第107页：分形图像压缩编码·程序代码

DD1(x,y)=k; DD2(x,y)=l; end

if e<0.05 %给定误差值0.005,如果误差小于该误差则可以不

进行循环，直接进行与下一个R相配的D的寻找

break; end

end %第三层循环结束

if e<0.05 break; end

end %第二层循环结束

if e<0.05 break; end end

end %第一层循环结束

end save lenadata DD1 DD2 nn r h %保存压缩系数

第108页：分形图像解码程序代码

2)分形图像解码

%load impdata; load lenadata; N=256; J=4; D=2*J; I=rand(N,N); I1=ones(N,N); times=10; for t=1:times for i=1:J:(N-J+1) for j=1:J:(N-J+1) x=floor((i-1)/J)+1; y=floor((j-1)/J)+1; B=out(I,D,DD1(x,y),DD2(x,y)); B1=average(B,D);

第109页：分形图像压缩编码·程序代码

R1=Ln(B1,J,nn(x,y)); R1=r(x,y)*R1; R2=out(I1,J,DD1(x,y),DD2(x,y)); R2=h(x,y)*R2; R=R1+R2; I=in(I,R,J,i,j); end end end imshow(I); save decodelena I; %将解码后的图像存储于dat文件中3) Out函数 function B=out(I,j,k,l);

%取块操作。该函数作用于图像I，从I上取出一个二维子块，子块大小为2(J)*2(J)；

%k，l表示子块在二维图像中的坐标，且在子块的左上角

B(1:j,1:j)=I(k:k+j-1,l:l+j-1);

第110页：分形图像压缩编码·程序代码

4) Average函数 function R=average(I,j);

%求平均，表示将2(J)*2(J)大小子块的灰度图像中的邻近四像素点灰度求平均抽取，变成

2(J-1)*2(J-1)大小的子块 for k=1:j/2 for l=1:j/2 R(k,l)=round((I(2*k-1,2*l-1)+I(2*k,2*l)+I(2*k-1,2*l)+I(2*k,2*l-1))/4); end end5) 8种仿射变换函数Ln function D=Ln(R,j,i);

%仿射变换，共8种仿射变换

for k=0:j-1 for l=0:j-1 if i==1 D(k+1,l+1)=R(k+1,l+1); elseif i==2 D(k+1,l+1)=R(k+1,j-l); elseif i==3 D(k+1,l+1)=R(j-k,l+1);

第111页：分形图像压缩编码·程序代码

elseif i==4 D(k+1,l+1)=R(l+1,k+1); elseif i==5 D(k+1,l+1)=R(j-l,j-k); elseif i==6 D(k+1,l+1)=R(l+1,j-k); elseif i==7 D(k+1,l+1)=R(j-k,j-l); elseif i==8 D(k+1,l+1)=R(j-l,k+1); else D(k+1,l+1)=R(k+1,l+1); end end end

第112页：分形图像压缩编码·程序代码

6) Rhdp函数，比较两个块中元素灰度值的差别

function [r,h,dp]=rhdp(R,B,J);

%比较R、B中元素(灰度)值的差别，J为R、B矩阵的行数

m=J*J; A1=0; A2=0; A3=0; A4=0; A5=0; temp=ones(J); for i=1:J for j=1:J A1=A1+R(i,j)*B(i,j); A2=A2+R(i,j); A3=A3+B(i,j); A4=A4+B(i,j)*B(i,j); A5=A5+R(i,j)*R(i,j); end end if (R==temp)&(B==temp)

第113页：分形图像压缩编码·程序代码

r=1; h=0; dp=0; else if (m*A4-A3*A3)==0 %if (m*A4-A3*A3)==0,we do't select itas the last result.??? r=1; h=0; dp=100; else r=(m*A1-A2*A3)/(m*A4-A3*A3); h=(A2-r*A3)/m; dp=(A5+r*(r*A4-2*A1+2*h*A3)+h*(h*m-2*A2))/m; end end

第114页：放块操作程序代码

7)放块操作

function M=in(I,B,j,k,l);

%放块操作，将2(J)+2(J)的二维子块插入全0的图中，(k，l)为子块左上角坐标

I(k:k+j-1,l:l+j-1)=B(1:j,1:j); M=I;

8)计算信噪比的函数

N=258; %the size of image block. J=4 lena=imread('lena1.tif'); old=double(lena); load decodelena; %decodedata4; new=I; ysb=J*J/5 S=0; No=0; for i=1:N for j=1:N S=S+old(i,j)^2; No=No+(new(i,j)-old(i,j))^2; end end SNR=10*log10(S/No)

第115页：分形图像压缩编码·运行结果与分析

5．运行结果与分析用以上所述方法，将大小为256×256×8的lena图像进行压缩及解压缩，结果如图5.39和图5.40所示。图5.39 lena原图图5.40解码后的lena图像

第116页：分形图像压缩编码·思考题

6．思考题

(1)分形图像压缩编码压缩后需要存储的内容是什么？

(2)为什么分形图像压缩编码的压缩时间很长？应该怎样改进？本章小结本章通过图像配准、图像去噪、直方图均衡、基于小波的图像压缩、基于分形理论的图像压缩等实例简介了MATLAB在图像处理中的应用。通过本章的学习，对于图像的读取、显示、保存等方法有一定的了解，并能学会在MATLAB中进行图像处理的基本方法。图像处理工具箱是MATLAB中非常强大的工具箱之一。MATLAB对变量以矩阵形式处理的方式使其非常适合图像处理中的应用。数字图像处理理论和知识非常广泛，本章只涉及其中很小的一部分。希望能起到抛砖引玉的作用。

第117页：结束页

Thanks书山有路勤为径，学海无涯苦作舟。好好学习，天天向上