数字图像处理实验报告-基于

1、华东师范大学电子工程系实验 1：图像灰度级修正【 实验目的 】掌握常用的图像灰度级修正方法(灰度变换法和直方图均衡化)，加深对直方图的理解。观察图像的增强效果，对灰度级修正前后的图像加以比较。【 实验内容 】1)编程实现图像的灰度变换，改变图像的输入、输出映射参数范围(线性拉伸和反比) ；2)修改参数gammOS (大于、小于、等于1),观察处理结果；3)对图像直方图作均衡化处理，显示均衡前后的图像及其直方图。【实验代码】original=imread( '' );linstr=imadjust(original, ,0 1);%线性拉伸opposite=imadjust(or

2、iginal,0 1,1 0);%反比above=imadjust(original,0 1,0 1,2);%gamma>1equal=imadjust(original,0 1,0 1,1);%gamma=1below=imadjust(original,0 1,0 1,;%gamma<1%对图像均衡化' 原图像 ' );' 线性拉伸' );' 反比 ' );gamma>1' );gamma=1' );gamma<1' );' 原图像直方图' );均衡后的图像' );均衡图

3、像的直方图' )subplot(3,3,1);imshow(original);title( subplot(3,3,2);imshow(linstr);title( subplot(3,3,3);imshow(opposite);title( subplot(3,3,4);imshow(above);title( subplot(3,3,5);imshow(equal);title( subplot(3,3,6);imshow(below);title( subplot(3,3,7);imhist(original);title( histequal=histeq(original)

4、;subplot(3,3,8);imshow(histequal);title( subplot(3,3,9);imhist(histequal);title( axis(0 256 0 2000);【输出图像】【实验思考】根据以下图片以及实验结果可知gamma>1寸图像整体变暗，灰度级整体变小;gammav时图像整体变亮，灰度级整体变小；而 gamma二时，图像维持不变。实验2：图像的平滑滤波【实验目的】平滑的目的是减少噪声对图像的影响。掌握线性滤波和中值滤波两种最典型、最常用的图像平滑方法，对输出结果加以比较、加深理解。【实验内容】1)编写并调试窗口为3X3、5X5的平滑滤波函数；如

5、1 1 1;1 1 1 ;1 1 1/9、 1 2 1;2 4 2;1 2 1/16 等)2)编写并调试窗口为3X 3、5X5的中值滤波函数。3)比较均值滤波和中值滤波的优缺点，分析窗口尺寸对滤波结果的影响。附：可供参考的Matlab 函数有 imnoise 、 imfilter 、 medfilt2【实验代码】function fliterI = imread( '' );%原始图像读取J = imnoise(I, 'salt & pepper' ,; %含噪图像加椒盐噪声subplot(2,3,1);imshow(J);title(' 含噪图

6、像' );Newbuf1=AverageFilter(J,256,256,3); subplot(2,3,2);imshow(Newbuf1);title(Newbuf2=AverageFilter(J,256,256,5); subplot(2,3,3);imshow(Newbuf2);title(W=1 2 1;2 4 2;1 2 1/16;Newbuf3=WeighFilter(J,W,256,256,3); subplot(2,3,4);imshow(Newbuf3);title(Newbuf4=MedianFilter(J,256,256,3); subplot(2,3,5)

7、;imshow(Newbuf4);title(Newbuf5=MedianFilter(J,256,256,5);subplot(2,3,6);imshow(Newbuf5);title(%3X 3标准平均，调用均值滤波函数 3 X 3标准平均');%5X 5标准平均，调用均值滤波函数 5 X 5标准平均');%设置加权平均掩膜%3X 3加权平均，调用加权平均函数3 X 3加权平均)；%3X 3中值滤波，调用中值滤波函数3 X 3中值滤波)；%5X 5中值滤波，调用中值滤波函数5 X 5中值滤波)；%标准平均滤波函数function Newbuf=AverageFilter(O

8、ldbuf,M,N,m)%Newbuf滤波后图像矩阵%Oldbuf含噪图像矩阵%M、 N含噪图像像素矩阵行、列%m均值滤波窗口大小f=zeros(M+m-1,N+m-1);%将原图像像素复制到f 矩阵上，空出(m-1)/2 大小的边界f(m-1)/2+1:M+(m-1)/2,(m-1)/2+1:N+(m-1)/2)=Oldbuf(1:M,1:N);%将与边界相邻的(m-1)/2 行(或列)的像素值复制到边界，以填充边界f(m-1)/2+1:M+(m-1)/2,1:(m-1)/2)=Oldbuf( : ,1:(m-1)/2);f(m-1)/2+1:M+(m-1)/2,N+(m-1)/2:N+m-

9、1)=Oldbuf( : ,N-(m-1)/2:N);f(1:(m-1)/2,(m-1)/2+1:N+(m-1)/2)=Oldbuf(1:(m-1)/2, : );f(M+(m-1)/2:M+m-1,(m-1)/2+1:N+(m-1)/2)=Oldbuf(M-(m-1)/2:M, : );g=zeros(M+m-1,N+m-1);Im=zeros(M,N);%根据公式计算出处理后g(x,y) 的像素值for x=(m-1)/2+1:M+(m-1)/2for y=(m-1)/2+1:N+(m-1)/2for s=-(m-1)/2:(m-1)/2for t=-(m-1)/2:(m-1)/2g(x,

10、y)=g(x,y)+f(x+s,y+t)*1/(m*m);endendendendIm(1:M,1:N)=g(m-1)/2+1:M+(m-1)/2,(m-1)/2+1:N+(m-1)/2);%t-double型转换为uint8型才可以用imshow正常显示 Newbuf=uint8(Im);%加权平均滤波函数function Newbuf=WeighFilter(Oldbuf,W,M,N,m)%Newbuf滤波后图像矩阵%Oldbuf含噪图像矩阵%W掩模%M、 N含噪图像像素矩阵行、列%m掩模模板窗口大小f=zeros(M+m-1,N+m-1);%等原图像像素复制到f矩阵上，空出(m-1)/2

11、大小的边界f(m-1)/2+1:M+(m-1)/2,(m-1)/2+1:N+(m-1)/2)=Oldbuf(1:M,1:N);%将与边界相邻的(m-1)/2 行(或列)的像素值复制到边界，以填充边界f(m-1)/2+1:M+(m-1)/2,1:(m-1)/2)=Oldbuf( : ,1:(m-1)/2);f(m-1)/2+1:M+(m-1)/2,N+(m-1)/2:N+m-1)=Oldbuf( : ,N-(m-1)/2:N);f(1:(m-1)/2,(m-1)/2+1:N+(m-1)/2)=Oldbuf(1:(m-1)/2, : );f(M+(m-1)/2:M+m-1,(m-1)/2+1:N+

12、(m-1)/2)=Oldbuf(M-(m-1)/2:M, : );g=zeros(M+m-1,N+m-1);Im=zeros(M,N);%根据公式计算出处理后g(x,y) 的像素值for x=(m-1)/2+1:M+(m-1)/2for y=(m-1)/2+1:N+(m-1)/2for s=-(m-1)/2:(m-1)/2for t=-(m-1)/2:(m-1)/2g(x,y)=g(x,y)+W(s+(m+1)/2,t+(m+1)/2)*f(x+s,y+t);endendendendIm(1:M,1:N)=g(m-1)/2+1:M+(m-1)/2,(m-1)/2+1:N+(m-1)/2);Ne

13、wbuf=uint8(Im);%中值滤波函数 function Newbuf=MedianFilter(Oldbuf,M,N,m)%Newbuf滤波后图像矩阵%Oldbuf含噪图像矩阵%M、 N含噪图像矩阵像素行、列%m中值滤波窗口大小f=zeros(M+m-1,N+m-1);%将原图像像素复制到f 矩阵上，空出(m-1)/2 大小的边界f(m-1)/2+1:M+(m-1)/2,(m-1)/2+1:N+(m-1)/2)=Oldbuf(1:M,1:N);%将与边界相邻的(m-1)/2 行(或列)的像素值复制到边界，以填充边界f(m-1)/2+1:M+(m-1)/2,1:(m-1)/2)=Oldb

14、uf( : ,1:(m-1)/2);f(m-1)/2+1:M+(m-1)/2,N+(m-1)/2:N+m-1)=Oldbuf( : ,N-(m-1)/2:N);f(1:(m-1)/2,(m-1)/2+1:N+(m-1)/2)=Oldbuf(1:(m-1)/2, : );f(M+(m-1)/2:M+m-1,(m-1)/2+1:N+(m-1)/2)=Oldbuf(M-(m-1)/2:M, : );g=zeros(M+m-1,N+m-1);Im=zeros(M,N);for x=(m-1)/2+1:M+(m-1)/2for y=(m-1)/2+1:N+(m-1)/2j=1;for s=-(m-1)/

15、2:(m-1)/2for t=-(m-1)/2:(m-1)/2a(j)=f(x+s,y+t);%将窗口里的二维元素变成一维元素j=j+1;endendg(x,y)=SeekMid(a,m);endendIm(1:M,1:N)=g(m-1)/2+1:M+(m-1)/2,(m-1)/2+1:N+(m-1)/2);Newbuf=uint8(Im);%找出中值函数function mid=SeekMid(winbuf,m)%mid排序后的中值%winbuf待排序窗口%m窗口大小%采用冒泡排序方法将窗口像素值从小到大排列，返回中间像素值for i=1:m*m-1for k=1:m*m-iif winbu

16、f(k)>winbuf(k+1)temp=winbuf(k);winbuf(k)=winbuf(k+1);winbuf(k+1)=temp;endendendmid=winbuf(m*m+1)/2);【输出图像】【实验思考】1. 比较均值滤波和中值滤波的优缺点均值滤波可以减小图像灰度级的“ sharp”变化，可以降低噪声，但是降噪的同时也使边缘部分变得模糊，还可以平滑伪轮廓，去除图像中的不相关的小于掩模尺寸的细节。中值滤波器的主要功能是使具有不同灰度的点看起来更接近它的相邻点，去除那些相对于其邻域像素更亮或更暗、且区域小于n2/2 的孤立像素集。中值滤波对降低某些类型的随机噪声性能优异，

17、模糊程度低。在处理椒盐噪声时，均值滤波使图像变得模糊，并且噪声去除性能很差，而中值滤波的效果却很好。显然，中值滤波比均值滤波更适合去除椒盐噪声。2. 分析窗口尺寸对滤波结果的影响窗口尺寸越大，图像越模糊，图像边缘和与掩膜大小接近的细节受到的影响也越大实验 3：图像的锐化处理【实验目的】锐化的目的是加强图像的边界和细节，熟悉 Robert、 Sobel 和 Laplace 算子进行检测，使图像特征(如边缘、轮廓等)进一步增强并突出。【实验内容】1)编写Robert 算子滤波函数；2)编写Sobel 算子滤波函数；3)编写Laplace算子滤波函数；4)编写限幅和标定函数，给出增强后的图像。【实验

18、代码】function EX3I=imread( '' );subplot(2,4,1);imshow(I);title(' 原始图像' );rob=RobertFilter(I);Robert 算子滤波结果' );Robert 算子增强结果' );subplot(2,4,2);imshow(rob);title(R1=I+rob;la1=LimitAmplitude(R1);subplot(2,4,6);imshow(la1);title( a2=-1 -2 -1;0 0 0;1 2 1;b2=-1 0 1;-2 0 2;-1 0 1;sob=

19、SobelFilter(I,a2,b2);subplot(2,4,3);imshow(sob);title(Sobel 算子滤波结果' );R2=I+sob;la2=LimitAmplitude(R2);subplot(2,4,7);imshow(la2);title('Sobel 算子增强结果' );% s=0 1 0;1 -4 1;0 1 0;s=1 1 1;1 -8 1;1 1 1;lap=LapFilter(I,s);cal=Calibration(lap);subplot(2,4,4);imshow(cal);title(Laplace 算子滤波结果'

20、 );lap=uint8(lap);lapr=I-lap;lapr3=LimitAmplitude(lapr);subplot(2,4,8);imshow(lapr3);title(Laplace 算子增强结果' );%Robert算子滤波function rob=RobertFilter(F)a1=-1 0;0 1;b1=0 -1;1 0;%Robert 算子模板M,N=size(F);f=zeros(M+1,N+1);f(1:M,1:N)=F(1:M,1:N);f(1:M,N+1:N+1)=F( : ,N:N);f(M+1:M+1,1:N)=F(M:M, : );%边界填充g=ze

21、ros(M+1,N+1);for x=1:Mfor y=1:Nmod=f(x,y) f(x,y+1);f(x+1,y) f(x+1,y+1);gsx=a1.*mod;gsy=b1.*mod;g(x,y)=abs(sum(gsx(:)+abs(sum(gsy(:); endendIm=zeros(M,N);Im(1:M,1:N)=g(1:M,1:N);rob=uint8(Im);%Sobel算子滤波function sob=SobelFilter(F,sx,sy)%sx,sy 为 Sobel 算子模板 M,N=size(F);m,n=size(sx);f=zeros(M+m-1,N+n-1);f

22、(m-1)/2+1:M+(m-1)/2,(n-1)/2+1:N+(n-1)/2)=F(1:M,1:N);f(m-1)/2+1:M+(m-1)/2,1:(n-1)/2)=F( : ,1:(n-1)/2);f(m-1)/2+1:M+(m-1)/2,N+(n-1)/2:N+m-1)=F( : ,N-(n-1)/2:N);f(1:(m-1)/2,(n-1)/2+1:N+(n-1)/2)=F(1:(m-1)/2, : );f(M+(m-1)/2:M+m-1,(n-1)/2+1:N+(n-1)/2)=F(M-(m-1)/2:M, : );%边界填充g=zeros(M+m-1,N+n-1);for x=(m

23、-1)/2+1:M+(m-1)/2for y=(n-1)/2+1:N+(n-1)/2mod=f(x-1,y-1) f(x-1,y) f(x-1,y+1); f(x,y-1) f(x,y) f(x,y+1);f(x+1,y-1)f(x+1,y) f(x+1,y+1);gsx=sx.*mod;gsy=sy.*mod;g(x,y)=abs(sum(gsx(:)+abs(sum(gsy(:); endendIm=zeros(M,N);Im(1:M,1:N)=g(m-1)/2+1:M+(m-1)/2,(n-1)/2+1:N+(n-1)/2);sob=uint8(Im);%Laplace 算子滤波func

24、tion lap=LapFilter(F,S)M,N=size(F);m,n=size(S);f=zeros(M+m-1,N+n-1);f(m-1)/2+1:M+(m-1)/2,(n-1)/2+1:N+(n-1)/2)=F(1:M,1:N);f(m-1)/2+1:M+(m-1)/2,1:(n-1)/2)=F( : ,1:(n-1)/2);f(m-1)/2+1:M+(m-1)/2,N+(n-1)/2:N+m-1)=F( : ,N-(n-1)/2:N);f(1:(m-1)/2,(n-1)/2+1:N+(n-1)/2)=F(1:(m-1)/2, : );f(M+(m-1)/2:M+m-1,(n-1)

25、/2+1:N+(n-1)/2)=F(M-(m-1)/2:M, : );g=zeros(M+m-1,N+n-1);for x=(m-1)/2+1:M+(m-1)/2for y=(n-1)/2+1:N+(n-1)/2mod=f(x-1,y-1) f(x-1,y) f(x-1,y+1); f(x,y-1) f(x,y) f(x,y+1);f(x+1,y-1)f(x+1,y) f(x+1,y+1);gs=S.*mod;g(x,y)=sum(gs(:); endendIm=zeros(M,N);Im(1:M,1:N)=g(m-1)/2+1:M+(m-1)/2,(n-1)/2+1:N+(n-1)/2);l

26、ap=Im;%限幅函数function la=LimitAmplitude(F)f=uint8(F);M,N=size(f);for x=1:Mfor y=1:Nif f(x,y)>=255;f(x,y)=255;else if f(x,y)<=0f(x,y)=0;else f(x,y)=f(x,y);%将灰度值限定在0到 255之间endendendendla=f;%标定函数function cal=Calibration(F)F=double(F);M,N=size(F);m1=min(min(F);for x=1:Mfor y=1:Nfm(x,y)=F(x,y)-m1;end

27、endm2=max(max(fm);Fm=double(fm);for x=1:Mfor y=1:Nfs(x,y)=255*(Fm(x,y)/m2);endendcal=uint8(fs);【输出图像】实验4：图像的统计特性【实验目的】观察序列图像帧内、帧间差值信号的分布曲线，理解图像在空间域和频率域上的统计特性及其在压缩中的重要性。【实验内容】1)编写帧内统计函数，计算差值图像(同一行差值、同一列差值)计分布曲线；2)编写帧间统计函数，计算差值图像(相邻帧的差值)，观察统计分布曲线( cla0/1 或 girl0/1 ) 。附：可供参考的Matlab 函数有sum、cat、 plot【实验代

28、码】function EX4oldbuf=imread( '' );I1=imread( '' );I2=imread( '' );newbuf1=Intrah(oldbuf,1);%帧内水平差值统计特性newbuf2=Intrah(oldbuf,0);%帧内垂直差值统计特性newbuf3=Inter(I1,I2); 帧间统计特性' 原始图像' );' 水平差值统计特性' );' 垂直差值统计特性' );'CLA1' );'CLA2' );' 帧间统计特性&#

29、39; );subplot(2,3,1);imshow(oldbuf);title( subplot(2,3,2);draw(newbuf1);title( subplot(2,3,3);draw(newbuf2);title( subplot(2,3,4);imshow(I1);title( subplot(2,3,5);imshow(I2);title( subplot(2,3,6);draw(newbuf3);title(function newbuf=Intrah(oldbuf,pop2)%帧内统计函数oldbuf=double(oldbuf);M,N=size(oldbuf);%防止

30、溢出将数据类型从uint8 型转换为double 型newbuf=zeros(1,511); if pop2=1 for i=1:Mfor j=1:N-1dH=oldbuf(i,j)-oldbuf(i,j+1);% 帧内水平灰度差值newbuf(dH+256)=newbuf(dH+256)+1;endendelsefor i=1:M-1for j=1:NdV=oldbuf(i,j)-oldbuf(i+1,j);newbuf(dV+256)=newbuf(dV+256)+1;end%帧间统计函数function newbuf=Inter(oldbuf,oldbuf1) oldbuf=double

31、(oldbuf);oldbuf1=double(oldbuf1);M,N=size(oldbuf);newbuf=zeros(1,511);for i=1:Mfor j=1:Ndt=oldbuf(i,j)-oldbuf1(i,j);% 计算帧间差值newbuf(dt+256)=newbuf(dt+256)+1; endendfunction draw(D)D=D/sum(D);x=-255:255;plot(x,D);axis(-100 100 0 );% 为了显示效果好缩小坐标轴范围【输出图像】实验6：方块编码【实验目的】掌握方块编码的基本方法及压缩性能。【实验内容】1)编程实现子块为nxn

32、的方块编码算法；2)分别取n=4和8的方块尺寸进行实验，计算重建图像的PSN序口压缩比。【实验代码】1. 主程序I=imread( '' );M,N=size(I);subplot(1,3,1);imshow(I);title(' 原图像 ' ); % 显示原图像I=double(I);newbuf1=BtcCode(I,4);PSNR1,Cr1=Analyze(I,newbuf1,M,N,4);subplot(1,3,2);imshow(uint8(newbuf1);title('4*4BTC 重建图像 ,PSNR=',num2str(PSNR

33、1), ' 压缩比 =' ,num2str(Cr1);newbuf2=BtcCode(I,8);PSNR2,Cr2=Analyze(I,newbuf2,M,N,8);'8*8BTC 重建图subplot(1,3,3);imshow(uint8(newbuf2);title(像 ,PSNR=' ,num2str(PSNR2), ' 压缩比 =' ,num2str(Cr2);2.function outbuf=BtcBlock(inbuf,n)%btc 方块编码算法函数%inbuf 子块数组%n 方块尺寸%对每个子块的图像数据分别计算xt 、 a0、

34、 a1 值，再用分辨率分量inbuf=double(inbuf);%(a0,a1) 替代方块原来的数据，最后放入方块图像数组中并返回该数组temp=0; %总的像素值temp0=0; %小于阀值的总像素temp1=0; %大于阀值的总像素q=0; %大于阀值的像素的个数 m=n*n; for i=1:n for j=1:n temp=temp+inbuf(i,j); endendxt=temp/m;%平均像素值即阀值for i=1:nfor j=1:nif inbuf(i,j)<xttemp0=temp0+inbuf(i,j);%得出小于阀值的总像素elsetemp1=temp1+inb

35、uf(i,j);%得出大于阀值的总像素q=q+1;%大于阀值的像素个数endendendif q=ma0=uint8(temp0/(m-q);%得出小于阀值的像素值endif q=0a1=uint8(temp1/q);%得出大于阀值的像素值endfor i=1:nfor j=1:nif inbuf(i,j)<xt outbuf(i,j)=a0;elseoutbuf(i,j)=a1;endnewbuf=BtcCode(oldbuf,n)%调用方块编码算法函数，输出编码后的图像M,N=size(oldbuf);row_num=M/n; %子块行数col_num=N/n; %子块列数row_s

36、tart=(0:row_num)*n+1;%子块起始行row_end=(1:row_num)*n;%子块终止行col_start=(0:col_num-1)*n+1;%子块起始列col_end=(1:row_num)*n;%子块终止列for i=1:row_numfor j=1:col_numf=oldbuf(row_start(i):row_end(i),col_start(j):col_end(j);%此式太长为方便书写定义foldbuf(row_start(i):row_end(i),col_start(j):col_end(j)=BtcBlock(f,n);%将原图像分成一个个子块，在

37、原图像里一个个对这些子块进行编码，编码后的结果保存原图像里endendnewbuf=oldbuf;%编码后的图像4.function PSNR,Cr=Analyze(I1,I2,M,N,n)%十算重建图像的PSN庠口压缩比m=n*n;mse=sum(sum(I1-I2).A2)/(M*N);PSNR=10*log10(255A2)/mse);Cr=8/(1+2*8/m);end【输出图像】实验7:JPEGE缩编码【实验目的】掌握nxn块的DCTS像变换及频谱特点。熟悉JPEG基本系统的图像编解码方法。【实验内容】1 )编程实现nxn块DC侬换的图像频谱显示，块 DCTS数按照 Zig-Zag扫

38、描并取部分进行图像重建，计算图像的均方根误差RMSE显示误差图 像和误差直方图。2)对8X8块的DCTS数，采用JPEGR认的量化矩阵进行量化和 反量化，计算原图像与重建图像之间的均方根误差RMSE并显示误差图像。【实验代码】1. 主程序F=imread( '' );subplot(231);imshow(F);title( '?- i ?'); %1示原图像F=double(F);F=F-128; %将原图像减小一半便于处理%十算原图像的8X8块的DC僚数，并转换为可视频谱图以便观察dctfre=DctBlock(F,8);subplot(232);imsho

39、w(log(abs(dctfre)*5+1),);title('8*8DCT 频谱显示' );% 表示将原图像的最大最小值之间的范围整体映射到0255之间，即做限幅DCTch=10;n=8;I,e,rmse1=ZigIDCT(F,dctfre,DCTch,n);取，num2str(DCTch),'个DCT(数时的重建图像);' 差值直方图，RMSE=',num2str(rmse1);subplot(233);imshow(uint8(I);title(subplot(234);imhist(uint8(abs(e);title( scale=4;newb

40、uf,err,rmse2=QuanIQuan(F,dctfre,n,scale);subplot(235);imshow(uint8(newbuf);title('scale 为 ' ,num2str(scale), ' 时的重建图像' );subplot(236);imshow(uint8(abs(err),);title(' 量化误差图像，RMSE=',num2str(rmse2); 2. function I,e,rmse1=ZigIDCT(oldbuf,dctfre,DCTch,n) %oldbuf: 原始图像数据 %dctfre:DCT

41、 系数矩阵 %DCTch每个分块中需要保留的DC琮数个数%n:分块的大小%e：原图像与保留部分DC僚数后的重建图像之间的误差矩阵 %按Zig-Zag扫描顺序，根据 DCTch>数，只保留64个 % DCT(数中的前DCTchl系数，对修改后的DC僚数用逆DC较换重建图彳象,得到DC倭%换的压缩图像。计算重建图像的均方根误差RMSE;显示误差图像和误差直方图。zigzag = 1 2 6 7 15 16 28 293 5 8 14 17 27 30 434 9 13 18 26 31 42 4410 12 19 25 32 41 45 5411 20 24 33 40 46 53 5521

42、 23 34 39 47 52 56 6122 35 38 48 51 57 60 6236 37 49 50 58 59 63 64;烟置 z扫描顺序mask=zigzag<=DCTch; %艮据当前DCTchtt得到"Z"字扫描的逻辑值，mask为logic类型 %寸修改后的DC琮数用逆DC使换重建图像，得到 DC段换的压缩图像 D=dctmtx(n); I=blkproc(dctfre,n n,'P1*(x.*P2)*P3' ,D',maskbuf,D);%1为重建的压缩图像矩阵e=oldbuf-I;%e:原图像与保留部分DC琮数后的重建

43、图像之间的误差矩阵I=I+128;rmse1=RMSE(e); enddctfre = DctBlock(oldbuf,n)%块DC函数：根据给定的n值，计算原图像的nxn块的DC僚数，并转换为可视频谱图以便观 察% oldbuf原始图像数据% n分块的大小% dctfre DCT 系数矩阵D=dctmtx(n);%D1返回Nix N的DC侵换矩阵，矩阵 A的DC直换可用Dx Ax D'来计算dctfre=blkproc(oldbuf,n,n,'P1*x*P2' ,D,D');%D'为 D 的转置end4.function newbuf,e,rmse2=

44、QuanIQuan(oldbuf,dctfre,n,scale)%量化和反量化函数：木K据给定的默认JPEGt化表，%寸每个n x n块的DCT(数进行量化和反量化，显示量化误差图像及其直方图%oldbuf: 原始图像数据%dctfre:DCT 系数矩阵%n:分块的大小%scale; 量化系数z= 16 11 10 16 24 40 51 6112 12 14 19 26 58 60 5514 13 16 24 40 57 69 5614 17 22 29 51 87 80 6218 22 37 56 68 109 103 7724 35 55 64 81 104 113 9249 64 78

45、 87 103 121 120 10172 92 95 98 112 100 103 99;僦认 JPEGt化表Qvalue=blkproc(dctfre,n n, 'round(x./P1)' ,scale*z);%量化IQvalue=blkproc(Qvalue,n n,'x.*P1' ,scale*z);%反量化%寸经过量化和反量化后的矩阵进行逆DCT5换得到重建图像矩阵D=dctmtx(n);newbuf=blkproc(IQvalue,n n,'P1*x*P2' ,D',D);e=newbuf-oldbuf;%物量化误差矩阵rm

46、se2=RMSE(e); %求均方根误差newbuf=newbuf+128;end5.function rmse=RMSE(oldbuf)%求均方根误差M,N=size(oldbuf);e=oldbuf.A2;rmse=sqrt(sum(e(:)/(M*N);end【输出图像】实验 8：运动估计【实验目的】熟悉运动估计的块匹配(BMA算法原理，编程实现全搜索算法(三步搜索或钻石搜索算法)，了解运动估计在混合编码器中的作用。【实验内容】1)编写全搜索算法函数，将运动矢量叠加到当前帧上并显示输出；2)显示输出预测帧、残差帧和重建图像，计算预测帧的PSNR。附：可供参考的Matlab 函数有 hol

47、d 、 quiver【实验代码】1. 主程序imgI=imread( '' ); %定义参考帧imgP=imread( '' ); %定义当前帧subplot(231);imshow(imgI);title(' 参考帧' );subplot(232);imshow(imgP);title(' 当前帧' );imgI=double(imgI);imgP=double(imgP);mbSize=16; %块尺寸为16*16p=7; %搜索窗口为(2p+1)*(2p+1)%基于块的全motionVect,EScomputations,b

48、lk_center,costs=ME_ES(imgP,imgI,mbSize,p);搜索算法imgMV(motionVect,imgP,blk_center);%画运动矢量图imgComp=motionComp(imgI,motionVect,mbSize);%根据运动矢量计算预测帧，并传输残差帧psnr=imgPSNR(imgP,imgComp);%计算峰值信噪比subplot(234);imshow(uint8(imgComp);title(' 预测帧 ,PSNR=',num2str(psnr);imgErr=imgP-imgComp;%残差帧cal=Calibration

49、(imgErr);%标定，显示更好效果subplot(235);imshow(cal);title(' 残差帧 ' );ChongJian=imgComp+imgErr;%根据运动矢量指明的位置及残差帧重建图像subplot(236);imshow(uint8(ChongJian);title(' 重建帧 ' );2.function motionVect,EScomputations,blk_center,costs = ME_ES(imgP, imgI, mbSize, p)%function:FS 算法：全搜索(Full Search/Exhaustive

50、 Search )%img:当前帧%imgI: 参考帧%mbSize:MB 寸%p:搜索窗口大小(2p+1) X ( 2p+1)%motionVect: 整像素精度MV%EScomputations: 搜索每个宏块所需的平均点数row,col=size(imgP);blk_center=zeros(2,row*col/(mbSizeA2);就义每个宏块中心点位置motionVect=zeros(2,row*col/(mbSizeA2);costs=ones(2*p+1,2*p+1)*65537;computations=0;%搜索的点数之和mbCount=1;for i = 1:mbSize:

51、row-mbSize+1for j = 1:mbSize:col-mbSize+1for m=-p:pfor n=-p:pref_blk_row=i+m;ref_blk_col=j+n;%当前帧起始行搜索范围，步长是块数%当前帧起始列搜索范围，步长是块数%定义每个宏块运动矢量%如果参考块的行列范围的任意一点在已经搜索过的宏块之外，则跳过，搜索%参考帧搜索框起始行%参考帧搜索框起始列点在已经搜索过的宏块之外，则跳过，搜索下一点if(ref_blk_row<1|ref_blk_row+mbSize-1>row|ref_blk_col<1|ref_blk_col+mbSize-1>col) continue ;end%否则计算该点SAD直costs(m+p+1,n+p+1)=costSAD(imgP(i:i+mbSize-1,j:j+mbSize-1),imgI(ref_blk_row:ref_blkrow+mbSize-1,ref_blk_col:ref_blk_col+mbSize-1);computations=computations+1; endend%记录中心点行坐标%记录中心点列坐标%找出有最小代价的块的下标%垂直运动矢量%水平运动矢量%重新赋值blk_center(1,mbCount) = i+ mbSize/2-1;blk_cent