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

华东师范大学电子工程系2017.6实验1：图像灰度级修正【实验目的】掌握常用的图像灰度级修正方法（灰度变换法和直方图均衡化），加深对直方图的理解。观察图像的增强效果，对灰度级修正前后的图像加以比较。【实验内容】1）编程实现图像的灰度变换，改变图像的输入、输出映射参数范围（线性拉伸和反比）；2）修改参数gamma值（大于、小于、等于1），观察处理结果；3）对图像直方图作均衡化处理，显示均衡前后的图像及其直方图。【实验代码】original=imread(lena.bmp);linstr=imadjust(original,0.3 0.7,0 1); %线性拉伸opposite=imadjust(original,0 1,1 0); %反比above=imadjust(original,0 1,0 1,2); %gamma1equal=imadjust(original,0 1,0 1,1); %gamma=1below=imadjust(original,0 1,0 1,0.5); %gamma1);subplot(3,3,5);imshow(equal);title(gamma=1);subplot(3,3,6);imshow(below);title(gamma1时图像整体变暗，灰度级整体变小；gammawinbuf(k+1) temp=winbuf(k); winbuf(k)=winbuf(k+1); winbuf(k+1)=temp; end endendmid=winbuf(m*m+1)/2); 【输出图像】【实验思考】1. 比较均值滤波和中值滤波的优缺点 均值滤波可以减小图像灰度级的“sharp”变化，可以降低噪声，但是降噪的同时也使边缘部分变得模糊，还可以平滑伪轮廓，去除图像中的不相关的小于掩模尺寸的细节。中值滤波器的主要功能是使具有不同灰度的点看起来更接近它的相邻点，去除那些相对于其邻域像素更亮或更暗、且区域小于n2/2的孤立像素集。中值滤波对降低某些类型的随机噪声性能优异，模糊程度低。在处理椒盐噪声时，均值滤波使图像变得模糊，并且噪声去除性能很差，而中值滤波的效果却很好。显然，中值滤波比均值滤波更适合去除椒盐噪声。2. 分析窗口尺寸对滤波结果的影响窗口尺寸越大，图像越模糊，图像边缘和与掩膜大小接近的细节受到的影响也越大实验3：图像的锐化处理【实验目的】锐化的目的是加强图像的边界和细节，熟悉Robert、Sobel和Laplace算子进行检测，使图像特征（如边缘、轮廓等）进一步增强并突出。【实验内容】1）编写Robert算子滤波函数；2）编写Sobel算子滤波函数；3）编写Laplace算子滤波函数；4）编写限幅和标定函数，给出增强后的图像。【实验代码】function EX3I=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;0 0 0;1 2 1;b2=-1 0 1;-2 0 2;-1 0 1; 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=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算子滤波结果); 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=zeros(M+1,N+1);for x=1:M for y=1:N mod=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(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-1)/2+1:M+(m-1)/2 for y=(n-1)/2+1:N+(n-1)/2 mod=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算子滤波function 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)/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)/2 for y=(n-1)/2+1:N+(n-1)/2 mod=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; %限幅函数function la=LimitAmplitude(F)f=uint8(F);M,N=size(f);for x=1:M for y=1:N if f(x,y)=255; f(x,y)=255; else if f(x,y)=0 f(x,y)=0; else f(x,y)=f(x,y); %将灰度值限定在0到255之间 end end endendla=f; %标定函数function cal=Calibration(F)F=double(F);M,N=size(F);m1=min(min(F);for x=1:M for y=1:N fm(x,y)=F(x,y)-m1; endendm2=max(max(fm);Fm=double(fm);for x=1:M for y=1:N fs(x,y)=255*(Fm(x,y)/m2); endendcal=uint8(fs);【输出图像】实验4：图像的统计特性【实验目的】观察序列图像帧内、帧间差值信号的分布曲线，理解图像在空间域和频率域上的统计特性及其在压缩中的重要性。【实验内容】1）编写帧内统计函数，计算差值图像（同一行差值、同一列差值），观察统计分布曲线；2）编写帧间统计函数，计算差值图像（相邻帧的差值），观察统计分布曲线（cla0/1或girl0/1）。附：可供参考的Matlab函数有sum、cat、plot【实验代码】function EX4oldbuf=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(帧间统计特性); function newbuf=Intrah(oldbuf,pop2) %帧内统计函数oldbuf=double(oldbuf);M,N=size(oldbuf); %防止溢出将数据类型从uint8型转换为double型newbuf=zeros(1,511);if pop2=1 for i=1:M for j=1:N-1 dH=oldbuf(i,j)-oldbuf(i,j+1);% 帧内水平灰度差值 newbuf(dH+256)=newbuf(dH+256)+1; end endelse for i=1:M-1 for j=1:N dV=oldbuf(i,j)-oldbuf(i+1,j); newbuf(dV+256)=newbuf(dV+256)+1; end endend function newbuf=Inter(oldbuf,oldbuf1) %帧间统计函数oldbuf=double(oldbuf);oldbuf1=double(oldbuf1);M,N=size(oldbuf); newbuf=zeros(1,511);for i=1:M for j=1:N dt=oldbuf(i,j)-oldbuf1(i,j);% 计算帧间差值 newbuf(dt+256)=newbuf(dt+256)+1; endend function draw(D)D=D/sum(D);x=-255:255;plot(x,D);axis(-100 100 0 0.5);% 为了显示效果好缩小坐标轴范围【输出图像】实验6：方块编码【实验目的】掌握方块编码的基本方法及压缩性能。【实验内容】1）编程实现子块为nn的方块编码算法；2）分别取n4和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.mfunction outbuf=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;for i=1:n for j=1:n temp=temp+inbuf(i,j); endendxt=temp/m; %平均像素值即阀值for i=1:n for j=1:n if inbuf(i,j)xt temp0=temp0+inbuf(i,j); %得出小于阀值的总像素 else temp1=temp1+inbuf(i,j); %得出大于阀值的总像素 q=q+1; %大于阀值的像素个数 end endendif q=m a0=uint8(temp0/(m-q); %得出小于阀值的像素值endif q=0 a1=uint8(temp1/q); %得出大于阀值的像素值endfor i=1:n for j=1:n if inbuf(i,j)xt outbuf(i,j)=a0; else outbuf(i,j)=a1; end endend3.BtcCode.mfunction newbuf=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; %子块终止列for i=1:row_num for j=1:col_num f=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; %编码后的图像4. 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(2552)/mse);Cr=8/(1+2*8/m);end【输出图像】实验7：JPEG压缩编码【实验目的】掌握nn块的DCT图像变换及频谱特点。熟悉JPEG基本系统的图像编解码方法。【实验内容】 1）编程实现nn块DCT变换的图像频谱显示，块DCT系数按照Zig-Zag扫描并取部分进行图像重建，计算图像的均方根误差RMSE，显示误差图像和误差直方图。 2）对88块的DCT系数，采用JPEG默认的量化矩阵进行量化和反量化，计算原图像与重建图像之间的均方根误差RMSE、并显示误差图像。【实验代码】1.主程序F=imread(lena.bmp);subplot(231);imshow(F);title(); %显示原图像F=double(F);F=F-128; %将原图像减小一半便于处理%计算原图像的88块的DCT系数，并转换为可视频谱图以便观察dctfre=DctBlock(F,8);subplot(232);imshow(log(abs(dctfre)*5+1),);title(8*8DCT频谱显示);%表示将原图像的最大最小值之间的范围整体映射到0255之间，即做限幅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);2. 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 7 15 16 28 29 3 5 8 14 17 27 30 43 4 9 13 18 26 31 42 44 10 12 19 25 32 41 45 54 11 20 24 33 40 46 53 55 21 23 34 39 47 52 56 61 22 35 38 48 51 57 60 62 36 37 49 50 58 59 63 64; %设置z扫描顺序mask=zigzag=DCTch; %根据当前DCTch值得到“Z”字扫描的逻辑值,mask为logic类型%对修改后的DCT系数用逆DCT变换重建图像，得到DCT变换的压缩图像D=dctmtx(n);I=blkproc(dctfre,n n,P1*(x.*P2)*P3,D,maskbuf,D); %I为重建的压缩图像矩阵e=oldbuf-I; %e:原图像与保留部分DCT系数后的重建图像之间的误差矩阵I=I+128;rmse1=RMSE(e);end3.DctBlock.mfunction dctfre = DctBlock(oldbuf,n) %分块DCT函数：根据给定的n值，计算原图像的nn块的DCT系数，并转换为可视频谱图以便观察 % oldbuf 原始图像数据 % n 分块的大小 % dctfre DCT系数矩阵D=dctmtx(n); %D是返回NN的DCT变换矩阵，矩阵A的DCT变换可用DAD来计算dctfre=blkproc(oldbuf,n,n,P1*x*P2,D,D); %D为D的转置end4. QuanIQuan.mfunction newbuf,e,rmse2=QuanIQuan(oldbuf,dctfre,n,scale)%量化和反量化函数：根据给定的默认JPEG量化表，%对每个nn块的DCT系数进行量化和反量化，显示量化误差图像及其直方图。%oldbuf:原始图像数据%dctfre:DCT系数矩阵%n:分块的大小%scale;量化系数z= 16 11 10 16 24 40 51 61 12 12 14 19 26 58 60 55 14 13 16 24 40 57 69 56 14 17 22 29 51 87 80 62 18 22 37 56 68 109 103 77 24 35 55 64 81 104 113 92 49 64 78 87 103 121 120 101 72 92 95 98 112 100 103 99; %默认JPEG量化表 Qvalue=blkproc(dctfre,n n,round(x./P1),scale*z); %量化 IQvalue=blkproc(Qvalue,n n,x.*P1,scale*z); %反量化 %对经过量化和反量化后的矩阵进行逆DCT变换得到重建图像矩阵 D=dctmtx(n); newbuf=blkproc(IQvalue,n n,P1*x*P2,D,D); e=newbuf-oldbuf; %e为量化误差矩阵 rmse2=RMSE(e); %求均方根误差 newbuf=newbuf+128;end5. RMSE.mfunction rmse=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算法：全搜索（Full Search/Exhaustive Search）%img:当前帧%imgI:参考帧%mbSize:MB尺寸%p:搜索窗口大小（2p+1）（2p+1）%motionVect:整像素精度MV%EScomputations:搜索每个宏块所需的平均点数row,col=size(imgP);blk_center=zeros(2,row*col/(mbSize2); %定义每个宏块中心点位置motionVect=zeros(2,row*col/(mbSize2); %定义每个宏块运动矢量costs=ones(2*p+1,2*p+1)*65537;computations=0; %搜索的点数之和mbCount=1;for i = 1:mbSize:row-mbSize+1 %当前帧起始行搜索范围，步长是块数 for j = 1:mbSize:col-mbSize+1 %当前帧起始列搜索范围，步长是块数 for m=-p:p for n=-p:p ref_blk_row=i+m; %参考帧搜索框起始行 ref_blk_col=j+n; %参考帧搜索框起始列 %如果参考块的行列范围的任意一点在已经搜索过的宏块之外，则跳过，搜索下一点 if (ref_blk_rowrow|ref_blk_colcol) 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; end end blk_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; %垂直运动矢量 motionVect(2,mbCount)=dy-p-1; %水平运动矢量 mbCount=mbCount+1; costs=ones(2*p+1,2*p+1)*65537; %重新赋值 endendEScomputations=computations/(mbCount-1); end3. motionComp.mfun