数字图像处理实验报告-基于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); %gamma>1

equal=imadjust(original,[0 1],[0 1],1); %gamma=1

below=imadjust(original,[0 1],[0 1],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([0 256 0 2000]);

【输出图像】

【实验思考】

根据以下图片以及实验结果可知gamma>1时图像整体变暗，

灰度级整体变小；gamma<1时图像整体变亮，灰度级整体变

小；而gamma=1时，图像维持不变。

实验2：图像的平滑滤波

【实验目的】

平滑的目的是减少噪声对图像的影响。掌握线性滤波和中值滤波两种最典型、最常用的图像平滑方法，对输出结果加以比较、加深理解。

【实验内容】

1）编写并调试窗口为3×3、5×5的平滑滤波函数；如[1 1 1;1 1 1 ;1 1 1]/9、[1 2 1;2 4 2;1 2 1]/16等）

2）编写并调试窗口为3×3、5×5的中值滤波函数。

3）比较均值滤波和中值滤波的优缺点，分析窗口尺寸对滤波结果的影响。附：可供参考的Matlab函数有imnoise、imfilter、medfilt2

【实验代码】

function fliter

I = imread('lena.bmp'); %原始图像读取

J = 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=[1 2 1;2 4 2;1 2 1]/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中值滤波');

%标准平均滤波函数

function Newbuf=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)的像素值

for x=(m-1)/2+1:M+(m-1)/2

for y=(m-1)/2+1:N+(m-1)/2

for s=-(m-1)/2:(m-1)/2

for t=-(m-1)/2:(m-1)/2

g(x,y)=g(x,y)+f(x+s,y+t)*1/(m*m);

end

end

end

end

Im(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);

%加权平均滤波函数

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大小的边界

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)的像素值

for x=(m-1)/2+1:M+(m-1)/2

for y=(m-1)/2+1:N+(m-1)/2

for s=-(m-1)/2:(m-1)/2

for t=-(m-1)/2:(m-1)/2

g(x,y)=g(x,y)+W(s+(m+1)/2,t+(m+1)/2)*f(x+s,y+t); end

end

end

end

Im(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 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)=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);

for x=(m-1)/2+1:M+(m-1)/2

for y=(m-1)/2+1:N+(m-1)/2

j=1;

for s=-(m-1)/2:(m-1)/2

for t=-(m-1)/2:(m-1)/2

a(j)=f(x+s,y+t);%将窗口里的二维元素变成一维元素

j=j+1;

end

end

g(x,y)=SeekMid(a,m);

end

end

Im(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-1

for k=1:m*m-i

if winbuf(k)>winbuf(k+1)

temp=winbuf(k);

winbuf(k)=winbuf(k+1);

winbuf(k+1)=temp;

end

end

end

mid=winbuf((m*m+1)/2);

【输出图像】

【实验思考】

1.比较均值滤波和中值滤波的优缺点

均值滤波可以减小图像灰度级的“sharp”变化，可以降低噪声，但是降噪的同时也使边缘部分变得模糊，还可以平滑伪轮廓，去除图像中的不相关的小于掩模尺寸的细节。

中值滤波器的主要功能是使具有不同灰度的点看起来更接近它的相邻点，去除那些相对于其邻域像素更亮或更暗、且区域小于n2/2的孤立像素集。中值滤波对降低某些类型的随机噪声性能优异，模糊程度低。

在处理椒盐噪声时，均值滤波使图像变得模糊，并且噪声去除性能很差，而中值滤波的效果却很好。显然，中值滤波比均值滤波更适合去除椒盐噪声。

2. 分析窗口尺寸对滤波结果的影响

窗口尺寸越大，图像越模糊，图像边缘和与掩膜大小接近的细节受到的影响也越大

实验3：图像的锐化处理

【实验目的】

锐化的目的是加强图像的边界和细节，熟悉Robert、Sobel和Laplace算子进行检测，使图像特征（如边缘、轮廓等）进一步增强并突出。

【实验内容】

1）编写Robert算子滤波函数；

2）编写Sobel算子滤波函数；

3）编写Laplace算子滤波函数；

4）编写限幅和标定函数，给出增强后的图像。

【实验代码】

function EX3

I=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(:)));

end

end

Im=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(:)));

end

end

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算子滤波

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(:));

end

end

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);

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;

elseif f(x,y)<=0

f(x,y)=0;

else f(x,y)=f(x,y);%将灰度值限定在0到255之间

end

end

end

end

la=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;

end

end

m2=max(max(fm));

Fm=double(fm);

for x=1:M

for y=1:N

fs(x,y)=255*(Fm(x,y)/m2);

end

end

cal=uint8(fs);

【输出图像】

实验4：图像的统计特性

【实验目的】

观察序列图像帧内、帧间差值信号的分布曲线，理解图像在空间域和频率域上的统计特性及其在压缩中的重要性。

【实验内容】

1）编写帧内统计函数，计算差值图像（同一行差值、同一列差值），观察统计分布曲线；

2）编写帧间统计函数，计算差值图像（相邻帧的差值），观察统计分布曲线（cla0/1

或girl0/1）。

附：可供参考的Matlab函数有sum、cat、plot

【实验代码】

function EX4

oldbuf=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

end

else

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

end

end

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;

end

end

function draw(D)

D=D/sum(D);

x=-255:255;

plot(x,D);

axis([-100 100 0 0.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.m

function 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);

end

end

xt=temp/m;%平均像素值即阀值

for i=1:n

for j=1:n

if inbuf(i,j)

temp0=temp0+inbuf(i,j);%得出小于阀值的总像素

else

temp1=temp1+inbuf(i,j);%得出大于阀值的总像素

q=q+1;%大于阀值的像素个数

end

end

end

if q~=m

a0=uint8(temp0/(m-q));%得出小于阀值的像素值

end

if q~=0

a1=uint8(temp1/q); %得出大于阀值的像素值

end

for i=1:n

for j=1:n

if inbuf(i,j)

outbuf(i,j)=a0;

else

outbuf(i,j)=a1;

end

end

end

3.BtcCode.m

function 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) ;%将原图像分成一个个子块，在原图像里一个个对这些子块进行编码，编码后的结果保存原图像里

end

end

newbuf=oldbuf;%编码后的图像

4. Analyze.m

function [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)]);

2. ZigIDCT.m

function [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 1

4 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);

end

3.DctBlock.m

function dctfre = 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的转置

end

4. QuanIQuan.m

function [newbuf,e,rmse2]=QuanIQuan(oldbuf,dctfre,n,scale)

%量化和反量化函数：根据给定的默认JPEG量化表，

%对每个n×n块的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;

end

5. RMSE.m

function 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*16

p=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.m

function [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/(mbSize^2));%定义每个宏块中心点位置

motionVect=zeros(2,row*col/(mbSize^2));%定义每个宏块运动矢量

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_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_b lk_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;%重新赋值

end

end

EScomputations=computations/(mbCount-1);

end

3. motionComp.m

function [imgComp] = motionComp(imgI, motionVect, mbSize)

%求预测帧的函数：由给定的运动矢量进行运动补偿计算预测帧

%imgI:参考帧

%motionVect:MV（dx为垂直分量，dy为水平分量）

%mbSize:MB尺寸

%imgComp:运动补偿后的图像

[row,col]=size(imgI);

mbCount=1;

for i = 1:mbSize:row-mbSize+1

for j = 1:mbSize:col-mbSize+1

ref_blk_row=i+motionVect(1,mbCount);%参考帧搜索块起始行

ref_blk_col=j+motionVect(2,mbCount);%参考帧搜索块起始列

imgComp(i:i+mbSize-1,j:j+mbSize-1)=imgI(ref_blk_row:ref_blk_row+mbSiz e-1,ref_blk_col:ref_blk_col+mbSize-1);

mbCount=mbCount+1;

end

end

end

4. minCost.m

function [minc,dx,dy] = minCost(costs)

%求具有最小SAD值的函数：找出具有最小代价的块的下标

%costs:包含当前宏块所有运动估计误差代价的SAD矩阵

%dx:MV的垂直分量

%dy:MV的水平分量

[row,col]=size(costs);

minc=costs(1,1);

dx=1;

dy=1;

for i=1:row

for j=1:col

if(costs(i,j)

minc=costs(i,j);

dx=i;

dy=j;

end

end

end

end

5. costSAD.m

function cost = costSAD(currentBlk,refBlk)

%function:SAD计算函数：对给定的两个块计算SAD值

%currentBlk：当前块

%refBlk：参考块

%mbSize：MB尺寸

%cost：两个块之间的误差代价（SAD）

cost=sum(sum(abs(currentBlk-refBlk)));

end