数字图像管理实验报告-根据Matlab

.-华东师范大学电子工程系2017.6.-实验1：图像灰度级修正【实验目的】掌握常用的图像灰度级修正方法（灰度变换法和直方图均衡化），加深对直方图的理解。观察图像的增强效果，对灰度级修正前后的图像加以比较。【实验内容】精品文档放心下载1）编程实现图像的灰度变换，改变图像的输入、输出映射参数范围（线性拉伸和反比）；精品文档放心下载2）修改参数gamma值（大于、小于、等于1），观察处理结果；精品文档放心下载3）对图像直方图作均衡化处理，显示均衡前后的图像及其直方图。感谢阅读【实验代码】original=imread('lena.bmp');linstr=imadjust(original,[0.30.7],[01]);%线性拉伸opposite=imadjust(original,[01],[10]);%反比above=imadjust(original,[01],[01],2);%gamma>1equal=imadjust(original,[01],[01],1);%gamma=1below=imadjust(original,[01],[01],0.5);%gamma<1谢谢阅读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('gamma>1');subplot(3,3,5);imshow(equal);title('gamma=1');subplot(3,3,6);imshow(below);title('gamma<1');subplot(3,3,7);imhist(original);title('原图像直方图');谢谢阅读histequal=histeq(original);%对图像均衡化subplot(3,3,8);imshow(histequal);title('均衡后的图像');subplot(3,3,9);imhist(histequal);title('均衡图像的直方图');精品文档放心下载axis([025602000]);【输出图像】.-【实验思考】根据以下图片以及实验结果可知gamma>1时图像整体变暗，灰度级整体变小；gamma<1时图像整体变亮，灰度级整体变小；而gamma=1时，图像维持不变。感谢阅读实验2：图像的平滑滤波【实验目的】平滑的目的是减少噪声对图像的影响。掌握线性滤波和中值滤波两种最典型、最常用的图像平滑方法，对输出结果加以比较、加深理解。感谢阅读【实验内容】1）编写并调试窗口为3×3、5×5的平滑滤波函数；如[111;111;111]/9、[1感谢阅读21;242;121]/16等）2）编写并调试窗口为3×3、5×5的中值滤波函数。感谢阅读3）比较均值滤波和中值滤波的优缺点，分析窗口尺寸对滤波结果的影响。感谢阅读附：可供参考的Matlab函数有imnoise、imfilter、medfilt2精品文档放心下载【实验代码】functionfliter=imread('lena.bmp');%原始图像读取精品文档放心下载=imnoise(I,'salt&pepper',0.02);%含噪图像加椒盐噪声精品文档放心下载subplot(2,3,1);imshow(J);title('含噪图像');谢谢阅读Newbuf1=AverageFilter(J,256,256,3);%3×3标准平均，调用均值滤波函数subplot(2,3,2);imshow(Newbuf1);title('3×3标准平均');精品文档放心下载Newbuf2=AverageFilter(J,256,256,5);%5×5标准平均，调用均值滤波函数subplot(2,3,3);imshow(Newbuf2);title('5×5标准平均');谢谢阅读W=[121;242;121]/16; %设置加权平均掩膜谢谢阅读Newbuf3=WeighFilter(J,W,256,256,3);%3×3加权平均，调用加权平均函数subplot(2,3,4);imshow(Newbuf3);title('3×3加权平均');感谢阅读Newbuf4=MedianFilter(J,256,256,3);%3×3中值滤波，调用中值滤波函数subplot(2,3,5);imshow(Newbuf4);title('3×3中值滤波');感谢阅读Newbuf5=MedianFilter(J,256,256,5);%5×5中值滤波，调用中值滤波函数subplot(2,3,6);imshow(Newbuf5);title('5×5中值滤波');精品文档放心下载%标准平均滤波函数functionNewbuf=AverageFilter(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)=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+(m-1)/2)=Oldbuf(M-(m-1)/2:M,:);谢谢阅读g=zeros(M+m-1,N+m-1);Im=zeros(M,N);%根据公式计算出处理后g(x,y)的像素值forx=(m-1)/2+1:M+(m-1)/2fory=(m-1)/2+1:N+(m-1)/2fors=-(m-1)/2:(m-1)/2fort=-(m-1)/2:(m-1)/2g(x,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);%将double型转换为uint8型才可以用imshow正常显示谢谢阅读Newbuf=uint8(Im);%加权平均滤波函数functionNewbuf=WeighFilter(Oldbuf,W,M,N,m)精品文档放心下载%Newbuf 滤波后图像矩阵%Oldbuf 含噪图像矩阵%W 掩模%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-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)的像素值forx=(m-1)/2+1:M+(m-1)/2fory=(m-1)/2+1:N+(m-1)/2fors=-(m-1)/2:(m-1)/2fort=-(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);Newbuf=uint8(Im);感谢阅读%中值滤波函数functionNewbuf=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)=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+(m-1)/2)=Oldbuf(M-(m-1)/2:M,:);感谢阅读g=zeros(M+m-1,N+m-1);Im=zeros(M,N);forx=(m-1)/2+1:M+(m-1)/2fory=(m-1)/2+1:N+(m-1)/2j=1;fors=-(m-1)/2:(m-1)/2fort=-(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);.-%找出中值函数functionmid=SeekMid(winbuf,m)精品文档放心下载%mid%winbuf%m

排序后的中值待排序窗口窗口大小%采用冒泡排序方法将窗口像素值从小到大排列，返回中间像素值谢谢阅读fori=1:m*m-1fork=1:m*m-iifwinbuf(k)>winbuf(k+1)temp=winbuf(k);winbuf(k)=winbuf(k+1);winbuf(k+1)=temp;endendendmid=winbuf((m*m+1)/2);【输出图像】【实验思考】比较均值滤波和中值滤波的优缺点均值滤波可以减小图像灰度级的“sharp”变化，可以降低噪声，但是降噪的同时也使边缘部分变得模糊，还可以平滑伪轮廓，去除图像中的不相关的小于掩模尺寸的细节。精品文档放心下载中值滤波器的主要功能是使具有不同灰度的点看起来更接近它的相邻点，去除那些相对于其邻域像素更亮或更暗、且区域小于n2/2的孤立像素集。中值滤波对降低某些类型的随机噪声性能优异，模糊程度低。感谢阅读在处理椒盐噪声时，均值滤波使图像变得模糊，并且噪声去除性能很差，而中值滤波的效果却很好。显然，中值滤波比均值滤波更适合去除椒盐噪声。谢谢阅读分析窗口尺寸对滤波结果的影响.-窗口尺寸越大，图像越模糊，图像边缘和与掩膜大小接近的细节受到的影响也越大感谢阅读实验3：图像的锐化处理【实验目的】锐化的目的是加强图像的边界和细节，熟悉Robert、Sobel和Laplace算子进行检测，使图像特征（如边缘、轮廓等）进一步增强并突出。谢谢阅读【实验内容】1）编写Robert算子滤波函数；2）编写Sobel算子滤波函数；3）编写Laplace算子滤波函数；4）编写限幅和标定函数，给出增强后的图像。【实验代码】functionEX3I=imread('rice.bmp');subplot(2,4,1);imshow(I);title('原始图像');精品文档放心下载rob=RobertFilter(I);subplot(2,4,2);imshow(rob);title('Robert算子滤波结果');感谢阅读R1=I+rob;la1=LimitAmplitude(R1);subplot(2,4,6);imshow(la1);title('Robert算子增强结果');谢谢阅读a2=[-1-2-1;000;121];感谢阅读b2=[-101;-202;-101];精品文档放心下载sob=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=[010;1-41;010];s=[111;1-81;111];lap=LapFilter(I,s);cal=Calibration(lap);subplot(2,4,4);imshow(cal);title('Laplace算子滤波结果');谢谢阅读lap=uint8(lap);lapr=I-lap;lapr3=LimitAmplitude(lapr);谢谢阅读subplot(2,4,8);imshow(lapr3);title('Laplace算子增强结果');感谢阅读.-%Robert算子滤波functionrob=RobertFilter(F)精品文档放心下载a1=[-10;01];b1=[0-1;10];%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=zeros(M+1,N+1);forx=1:Mfory=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算子滤波functionsob=SobelFilter(F,sx,sy)感谢阅读%sx,sy为Sobel算子模板[M,N]=size(F);[m,n]=size(sx);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)/2+1:N+(n-1)/2)=F(M-(m-1)/2:M,:);%边界填充谢谢阅读g=zeros(M+m-1,N+n-1);forx=(m-1)/2+1:M+(m-1)/2fory=(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(:)));精品文档放心下载endend.-Im=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算子滤波functionlap=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)/2+1:N+(n-1)/2)=F(M-(m-1)/2:M,:);感谢阅读g=zeros(M+m-1,N+n-1);forx=(m-1)/2+1:M+(m-1)/2fory=(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);谢谢阅读lap=Im;%限幅函数functionla=LimitAmplitude(F)谢谢阅读f=uint8(F);[M,N]=size(f);forx=1:Mfory=1:Niff(x,y)>=255;f(x,y)=255;elseiff(x,y)<=0f(x,y)=0;elsef(x,y)=f(x,y);%将灰度值限定在0到255之间end精品文档放心下载endendendla=f;.-%标定函数functioncal=Calibration(F)谢谢阅读F=double(F);[M,N]=size(F);m1=min(min(F));forx=1:Mfory=1:Nfm(x,y)=F(x,y)-m1;endendm2=max(max(fm));Fm=double(fm);forx=1:Mfory=1:Nfs(x,y)=255*(Fm(x,y)/m2);endendcal=uint8(fs);【输出图像】实验4：图像的统计特性【实验目的】观察序列图像帧内、帧间差值信号的分布曲线，理解图像在空间域和频率域上的统计特性及其在压缩中的重要性。精品文档放心下载【实验内容】1）编写帧内统计函数，计算差值图像（同一行差值、同一列差值），观察统计精品文档放心下载分布曲线；2）编写帧间统计函数，计算差值图像（相邻帧的差值），观察统计分布曲线（cla0/1谢谢阅读.-girl0/1）。附：可供参考的Matlab函数有sum、cat、plot精品文档放心下载【实验代码】functionEX4oldbuf=imread('rice.bmp');精品文档放心下载I1=imread('CLA1.bmp');I2=imread('CLA2.bmp');newbuf1=Intrah(oldbuf,1);%帧内水平差值统计特性newbuf2=Intrah(oldbuf,0);%帧内垂直差值统计特性newbuf3=Inter(I1,I2);帧间统计特性精品文档放心下载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('CLA1');subplot(2,3,5);imshow(I2);title('CLA2');subplot(2,3,6);draw(newbuf3);title('帧间统计特性');谢谢阅读functionnewbuf=Intrah(oldbuf,pop2)%帧内统计函数oldbuf=double(oldbuf);精品文档放心下载[M,N]=size(oldbuf);%防止溢出将数据类型从uint8型转换为double型newbuf=zeros(1,511);感谢阅读ifpop2==1fori=1:Mforj=1:N-1dH=oldbuf(i,j)-oldbuf(i,j+1);%帧内水平灰度差值精品文档放心下载newbuf(dH+256)=newbuf(dH+256)+1;精品文档放心下载endendelsefori=1:M-1forj=1:NdV=oldbuf(i,j)-oldbuf(i+1,j);精品文档放心下载newbuf(dV+256)=newbuf(dV+256)+1;谢谢阅读endendendfunctionnewbuf=Inter(oldbuf,oldbuf1)%帧间统计函数谢谢阅读oldbuf=double(oldbuf);oldbuf1=double(oldbuf1);[M,N]=size(oldbuf);newbuf=zeros(1,511);.-fori=1:Mforj=1:Ndt=oldbuf(i,j)-oldbuf1(i,j);%计算帧间差值谢谢阅读newbuf(dt+256)=newbuf(dt+256)+1;感谢阅读endendfunctiondraw(D)D=D/sum(D);x=-255:255;plot(x,D);axis([-10010000.5]);%为了显示效果好缩小坐标轴范围感谢阅读【输出图像】实验6：方块编码【实验目的】掌握方块编码的基本方法及压缩性能。【实验内容】1）编程实现子块为n×n的方块编码算法；2）分别取n＝4和8的方块尺寸进行实验，计算重建图像的PSNR和压缩比。精品文档放心下载【实验代码】1.主程序I=imread('lena.bmp');[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(PSNR1),'压缩比=',num2str(Cr1)]);精品文档放心下载newbuf2=BtcCode(I,8);[PSNR2,Cr2]=Analyze(I,newbuf2,M,N,8);subplot(1,3,3);imshow(uint8(newbuf2));title(['8*8BTC重建图谢谢阅读像,PSNR=',num2str(PSNR2),'压缩比=',num2str(Cr2)]);感谢阅读2.BtcBlock.mfunctionoutbuf=BtcBlock(inbuf,n)谢谢阅读%btc方块编码算法函数%inbuf子块数组%n方块尺寸%对每个子块的图像数据分别计算xt、a0、a1值，再用分辨率分量%(a0,a1)替代方块原来的数据，最后放入方块图像数组中并返回该数组感谢阅读inbuf=double(inbuf);temp=0;%总的像素值temp0=0;%小于阀值的总像素temp1=0;%大于阀值的总像素q=0;%大于阀值的像素的个数m=n*n;fori=1:nforj=1:ntemp=temp+inbuf(i,j);endendxt=temp/m;%平均像素值即阀值fori=1:nforj=1:nifinbuf(i,j)<xttemp0=temp0+inbuf(i,j);%得出小于阀值的总像素精品文档放心下载elsetemp1=temp1+inbuf(i,j);%得出大于阀值的总像素精品文档放心下载q=q+1;%大于阀值的像素个数endendendifq~=ma0=uint8(temp0/(m-q));%得出小于阀值的像素值精品文档放心下载endifq~=0a1=uint8(temp1/q);%得出大于阀值的像素值精品文档放心下载endfori=1:n.-forj=1:nifinbuf(i,j)<xtoutbuf(i,j)=a0;elseoutbuf(i,j)=a1;endendend3.BtcCode.mfunctionnewbuf=BtcCode(oldbuf,n)感谢阅读%调用方块编码算法函数，输出编码后的图像[M,N]=size(oldbuf);row_num=M/n;%子块行数col_num=N/n;%子块列数row_start=(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;%子块终止列感谢阅读fori=1:row_numforj=1:col_numf=oldbuf(row_start(i):row_end(i),col_start(j):col_end(j));%此式太长为方便书写定义f感谢阅读oldbuf(row_start(i):row_end(i),col_start(j):col_end(j))=BtcBlock(f,n)谢谢阅读;%将原图像分成一个个子块，在原图像里一个个对这些子块进行编码，编码后的结果保存原图像里谢谢阅读endendnewbuf=oldbuf;%编码后的图像Analyze.mfunction[PSNR,Cr]=Analyze(I1,I2,M,N,n)感谢阅读%计算重建图像的PSNR和压缩比m=n*n;mse=sum(sum((I1-I2).^2))/(M*N);谢谢阅读PSNR=10*log10((255^2)/mse);精品文档放心下载Cr=8/(1+2*8/m);end【输出图像】.-实验7：JPEG压缩编码【实验目的】掌握n×n块的DCT图像变换及频谱特点。熟悉JPEG基本系统的图像编解码方法。【实验内容】1）编程实现n×n块DCT变换的图像频谱显示，块DCT系数按照Zig-Zag扫描并取部分进行图像重建，计算图像的均方根误差RMSE，显示误差图像和误差直方图。精品文档放心下载2）对8×8块的DCT系数，采用JPEG默认的量化矩阵进行量化和反量化，计算原图像与重建图像之间的均方根误差RMSE、并显示误差图像。精品文档放心下载【实验代码】1.主程序F=imread('lena.bmp');subplot(231);imshow(F);title('Ô-ͼÏñ');%显示原图像感谢阅读F=double(F);F=F-128;%将原图像减小一半便于处理%计算原图像的8×8块的DCT系数，并转换为可视频谱图以便观察精品文档放心下载dctfre=DctBlock(F,8);subplot(232);imshow(log(abs(dctfre)*5+1),[]);title('8*8DCT频谱显示');感谢阅读%[]表示将原图像的最大最小值之间的范围整体映射到0~255之间，即做限幅精品文档放心下载DCTch=10;n=8;[I,e,rmse1]=ZigIDCT(F,dctfre,DCTch,n);subplot(233);imshow(uint8(I));title(['取',num2str(DCTch),'个DCT系数时的精品文档放心下载重建图像']);subplot(234);imhist(uint8(abs(e)));title(['差值直方图，感谢阅读RMSE=',num2str(rmse1)]);scale=4;.-[newbuf,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)]);谢谢阅读ZigIDCT.mfunction[I,e,rmse1]=ZigIDCT(oldbuf,dctfre,DCTch,n)精品文档放心下载%oldbuf:原始图像数据%dctfre:DCT系数矩阵%DCTch:每个分块中需要保留的DCT系数个数%n:分块的大小%e:原图像与保留部分DCT系数后的重建图像之间的误差矩阵精品文档放心下载按Zig-Zag扫描顺序，根据DCTch参数，只保留64个谢谢阅读DCT系数中的前DCTch个系数，对修改后的DCT系数用逆DCT变换重建图像，得到DCT变感谢阅读换的压缩图像。计算重建图像的均方根误差RMSE；显示误差图像和误差直方图。感谢阅读zigzag=[1 2 6 715162829精品文档放心下载581417273043913182631424412192532414554202433404653552334394752566122353848515760623637495058596364];%设置z扫描顺序谢谢阅读mask=zigzag<=DCTch;%根据当前DCTch值得到“Z”字扫描的逻辑值,mask为logic类型精品文档放心下载%对修改后的DCT系数用逆DCT变换重建图像，得到DCT变换的压缩图像精品文档放心下载D=dctmtx(n);I=blkproc(dctfre,[nn],'P1*(x.*P2)*P3',D',maskbuf,D);%I为重建的压缩图像谢谢阅读矩阵e=oldbuf-I;%e:原图像与保留部分DCT系数后的重建图像之间的误差矩阵感谢阅读I=I+128;rmse1=RMSE(e);end3.DctBlock.mfunctiondctfre=DctBlock(oldbuf,n)谢谢阅读%分块DCT函数：根据给定的n值，计算原图像的n×n块的DCT系数，并转换为可视频谱图以便观察谢谢阅读%oldbuf 原始图像数据% n 分块的大小% dctfre DCT系数矩阵D=dctmtx(n);%D是返回N×N的DCT变换矩阵，矩阵A的DCT变换可用D×A×D’来计算dctfre=blkproc(oldbuf,[n,n],'P1*x*P2',D,D');%D'为D的转置感谢阅读end4.QuanIQuan.mfunction[newbuf,e,rmse2]=QuanIQuan(oldbuf,dctfre,n,scale)谢谢阅读.-%量化和反量化函数：根据给定的默认JPEG量化表，%对每个n×n块的DCT系数进行量化和反量化，显示量化误差图像及其直方图。感谢阅读%oldbuf:原始图像数据%dctfre:DCT系数矩阵%n:分块的大小%scale;量化系数z=[1611101624405161121214192658605514131624405769561417222951878062182237566810910377243555648110411392496478871031211201017292959811210010399];%默认JPEG量化表Qvalue=blkproc(dctfre,[nn],'round(x./P1)',scale*z);%量化IQvalue=blkproc(Qvalue,[nn],'x.*P1',scale*z);%反量化精品文档放心下载%对经过量化和反量化后的矩阵进行逆DCT变换得到重建图像矩阵谢谢阅读D=dctmtx(n);newbuf=blkproc(IQvalue,[nn],'P1*x*P2',D',D);感谢阅读e=newbuf-oldbuf;%e为量化误差矩阵感谢阅读rmse2=RMSE(e);%求均方根误差newbuf=newbuf+128;end5.RMSE.mfunctionrmse=RMSE(oldbuf)感谢阅读%求均方根误差[M,N]=size(oldbuf);e=oldbuf.^2;rmse=sqrt(sum(e(:))/(M*N));精品文档放心下载end【输出图像】.-实验8：运动估计【实验目的】熟悉运动估计的块匹配（BMA）算法原理，编程实现全搜索算法（三步搜索或钻石搜索算法），了解运动估计在混合编码器中的作用。谢谢阅读【实验内容】1）编写全搜索算法函数，将运动矢量叠加到当前帧上并显示输出；谢谢阅读2）显示输出预测帧、残差帧和重建图像，计算预测帧的PSNR。谢谢阅读附：可供参考的Matlab函数有hold、quiver感谢阅读【实验代码】1.主程序imgI=imread('CLA1.bmp');%定义参考帧精品文档放心下载imgP=imread('CLA2.bmp');%定义当前帧谢谢阅读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,blk_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(imgErr);%标定，显示更好效果谢谢阅读subplot(235);imshow(cal);title('残差帧');感谢阅读ChongJian=imgComp+imgErr;%根据运动矢量指明的位置及残差帧重建图像subplot(236);imshow(uint8(ChongJian));title('重建帧');精品文档放心下载2.ME_ES.mfunction[motionVect,EScomputations,blk_center,costs]=ME_ES(imgP,感谢阅读imgI,mbSize,p)%function:FS算法：全搜索（FullSearch/ExhaustiveSearch）谢谢阅读%img:当前帧%imgI:参考帧%mbSize:MB尺寸.-%p:搜索窗口大小（2p+1）×（2p+1）%motionVect:整像素精度MV%EScomputations:搜索每个宏块所需的平均点数谢谢阅读[row,col]=size(imgP);blk_center=zeros(2,row*col/(mbSize^2));%定义每个宏块中心点位置motionVect=zeros(2,row*col/(mbSize^2));%定义每个宏块运动矢量感谢阅读costs=ones(2*p+1,2*p+1)*65537;感谢阅读computations=0;%搜索的点数之和mbCount=1;fori=1:mbSize:row-mbSize+1%当前帧起始行搜索范围，步长是块数forj=1:mbSize:col-mbSize+1%当前帧起始列搜索范围，步长是块数感谢阅读form=-p:pforn=-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_blk_row+mbSize-1,ref_blk_col:ref_blk_col+mbSize-1));谢谢阅读computations=computations+1;精品文档放心下载endendblk_center(1,mbCount)=i+mbSize/2-1;%记录中心点行坐标谢谢阅读blk_center(2,mbCount)=j+mbSize/2-1;%记录中心点列坐标精品文档放心下载[minc,dx,dy]=minCost(costs);%找出有最小代价的块的下标谢谢阅读motionVect(1,mbCount)=dx-p-1;%垂直运动矢量精品文