数字图像处理实验答案2.doc

实验名称：图像的锐化处理

一、 实验目的
学习用锐化处理技术来加强图像的目标边界和图像细节。对图像进行梯度算子、Roberts算子、Sobel算子边缘检测处理和Laplace算子边缘增强处理，使图像的某些特征（如边缘、轮廓等）得以进一步的增强和突出。

二、 实验内容
（1） 编写梯度算子和Roberts算子滤波函数。
（2） 编写Sobel算子滤波函数。
（3） 编写Laplace算子边缘增强滤波函数。
（4） 观察频域中用高低通滤波技术对图像进行的平滑和锐化处理。

三、 实验方法及编程
用Matlab进行编程，主要是对照着实验书上面的算法流程图进行代码的编 写，主要编写代码如下所示：
function newbuf=RobFilter(oldbuf,M,N);
% *************************************************************************
% 函数名称：
% RobFilter()
% 说明：
% 'Robert梯度'滤波算法。
% *************************************************************************
for i=1:M-1
for j=1:N-1
newbuf(i,j)=abs(oldbuf(i,j)-oldbuf(i+1,j+1))+abs(oldbuf(i+1,j)-oldbuf(i,j+1));
end
end

%-------------------------------------------------------------------------

function newbuf=SobFilter(oldbuf,M,N);
% *************************************************************************
% 函数名称：
% SobFilter()
% 说明：
% 'Sobel'滤波算法。
% *************************************************************************
for i=2:M-1
for j=2:N-1
sx=oldbuf(i+1,j-1)+2*oldbuf(i+1,j)+oldbuf(i+1,j+1)-oldbuf(i-1,j-1)-2*oldbuf(i-1,j)-oldbuf(i-1,j+1);
sy=oldbuf(i-1,j+1)+2*oldbuf(i,j+1)+oldbuf(i+1,j+1)-oldbuf(i-1,j-1)-2*oldbuf(i,j-1)-oldbuf(i+1,j-1);
newbuf(i,j)=abs(sx)+abs(sy);
end
end

%-------------------------------------------------------------------------
function newbuf=LapFilter(oldbuf,M,N);
% *************************************************************************
% 函数名称：
% LapFilter()
% 说明：
% 'Laplace'滤波算法。
% ************************************************************************
for i=2:M-1
for j=2:N-1
newbuf(i,j)=5*oldbuf(i,j)-oldbuf(i-1,j)-oldbuf(i+1,j)-oldbuf(i,j-1)-oldbuf(i,j+1);
end
end

四、 实验结果及分析（原图像和处理后的图像比较及分析）

如上图所示即为实验结果。如上图所示Robert算子能够检测出原图像的边缘，并且显示出来；Sobel在Robert算子的基础上增强了边缘，所以图像上结果会显示比原图像边缘要更加粗一些；Laplace算子与其他边缘增强方法不同，即其边缘的增强程度和边缘的方向无关，所以图像的细节明显比原来更加突出。

实验名称：图像方块编码

三、 实验目的
通过编程实验，掌握方块编码的基本方法及压缩性能。

四、 实验内容
编程实现子块为nxn的方块编码基本算法，分别取n=2,4,8方块尺寸进行编解码实验，计算编码后的均方误差和压缩比。
实验图像可为任意图像。

五、 实验方法及编程
用Matlab进行编程，主要是对照着实验书上面的算法流程图进行代码的编 写，下面是主要的实验代码：
function newbuf=BtcCode(oldbuf,M,N,n)
% *************************************************************************
% 函数名称：
% BtcCode() 方块编码函数
% 参数：
% oldbuf 原图像数组
% M N 原图像尺寸
% n 方块尺寸
% 说明：
% 调用方块编码算法函数，输出编码后的图像
% *************************************************************************
newbuf=oldbuf;
rowblocks=M/n;
colblocks=N/n;
for i=1:rowblocks-1
row=i*n;
for j=1:colblocks-1
col=j*n;
for i=1:n
for j=1:n
inbuf(i,j)=oldbuf(i+row,j+col);
end
end
outbuf=BtcBlock(inbuf,n);
for i=1:n
for j=1:n
newbuf(i+row,j+col)=outbuf(i,j);
end
end
end
end
%-------------------------------------------------------------------------
function outbuf=BtcBlock(inbuf,n)
% *************************************************************************
% 函数名称：
% btcblock() 方块编码算法函数
% 参数：
% inbuf 方块数组
% n 方块尺寸
% outbuf 存放处理后的方块图像
% 说明：
% 把原图像分成n*n子块，对每个方块的图像数据分别计算xt,a0,a1值，再用分辨率分
% 量(a0,a1)替代方块原来的数据 最后放入方块图像数组中并返回该数组。
% *************************************************************************
temp=0; temp0=0; temp1=0; q=0;m=n*n;
inbuf=double(inbuf);
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)>=xt)
q=q+1;
temp1=temp1+inbuf(i,j);
else
temp0=temp0+inbuf(i,j);
end
end
end
if q~=m
a0=round(temp0/(m-q));
end
if q~=0
a1=round(temp1/q);
end
for i=1:n
for j=1:n
if (inbuf(i,j)<xt)
outbuf(i,j)=a0;
else
outbuf(i,j)=a1;
end
end
End

六、 实验结果及分析（原图像和处理后的图像比较及分析）

如上图所示，即为实验结果。实验将图像分为m=n*n的子图像块，然后进行方块编码；由实验结果可知，窗口愈大时，方块效应越明显，方块编码图像变得越模糊；窗口愈小时，编码图像接近于原图像。

实验名称：图像线性预测编码

七、 实验目的
通过编程设计，掌握帧内DPCM的编解码方法（预测，量化）及其压缩性能。

八、 实验内容
（1） 编制一维前值预测DPCM编解码程序，预测系数取（1,0,0,0）。
（2） 编制二维前值预测DPCM编解码程序，预测系数取（1/2,1/4,0,1/4)。
（3） 重建图像与原图像误差图像，用绝对误差表示如下：

式中的n为放大因子，以便观察误差的分布情况。
本实验采用15个量化分层的主观量化器，其量化电平分别取（0，+5，+10，+17，+28，+39，+52，+67）。
分别计算重建后图像的PSNR。
实验的图像为cla0或cla1。

三、 实验方法及编程
用Matlab进行编程，主要是对照着实验书上面的算法流程图进行代码的编 写，主要的代码如下所示：
function newbuf=Dpcm_code(oldbuf,M,N,dim);
% *************************************************************************
% 函数名称：
% Dpcm_Code() “线性预测编解码器”算法函数
% 参数：
% oldbuf 原图像数组
% M,N 原图像尺度
% dim 选择预测编码维数
% newbuf 存放处理后的图像二维数组
% 说明：
% 根据线性预测编解码算法调用各个子模块，对原图像进行1D/2D线性预测编码和解码，
% 解码后的恢复图像放在newbuf数组中。
% *************************************************************************
global newbuf; %定义全局变量
for i=1:M
for j=1:N
pre_val=Predict_Value(i,j,N,dim);
err=oldbuf(i,j)-pre_val;
quan_err=Quant_Value(err);
res_val=Restor_Value(quan_err,pre_val);
newbuf(i,j)=Clip_Value(res_val);
end
end

%-------------------------------------------------------------------------
function newbuf=Error_Code(M,N,dim);
% *************************************************************************
% 函数名称：
% Error_Code() “传输误码解码器”算法函数
% 参数：
% M,N 原图像尺度
% dim 选择预测编码维数
% newbuf 存放处理后的图像二维数组
% 说明：
% 模拟信道传输过程中产生的误码，观察传输误码经解码后对恢复图像的影响，解码
% 后的恢复图像放在newbuf数组中。
% *************************************************************************
global newbuf;
wrong=zeros(M,N);
wrong(100,100)=120;
for i=1:M
for j=1:N
pre_val=Predict_Value(i,j,N,dim);
err=wrong(i,j);
quan_err=Quant_Value(err);
res_val=Restor_Value(quan_err,pre_val);
newbuf(i,j)=Clip_Value(res_val);
end
end

%-------------------------------------------------------------------------
function Pvalue=Predict_Value(row,col,N,dim)
% *************************************************************************
% {This function is used to give predicted value as linear predictor.
% The prediction formula is as follows:
% 1_D DPCM: ^x[i,j]=128 if j=1
% x'[i,j-1] if j>1
% 2_D DPCM: ^x[i,j]=128 if i=1,j=1
% x'[ i,j-1] if i=1,j>1
% x'[i-1,j ] if i>1,j=1 or N
% 1/2x'[ i,j-1] if i>1,j>1
% 1/8x'[i-1,j-1] (Pirsch's predictor)
% 1/4x'[i-1,j ]
% 1/8x'[i-1,j+1]
% Dim : Dimension of prediction
% Row : vertical coordinate of current pixel to be predicted
% COL : horizontal coordinate of current pixel to be predicted}
% *************************************************************************
global newbuf;

switch dim %预测算法编程
case 1
if col==1
Pvalue=128;
else
Pvalue=newbuf(row,col-1);
end
case 2
if (row==1&&col==1)
Pvalue=128;
end
if (row==1&&col>1)
Pvalue=newbuf(row,col-1);
end
if (row>1&&col==1)||(row>1&&col==N)
Pvalue=newbuf(row-1,col);
end
if (row>1&&col>1&&col<N)
Pvalue=(1/2)*newbuf(row,col-1)+(1/4)*newbuf(row-1,col-1)+...
(1/4)*newbuf(row-1,col+1);
end
end
%-------------------------------------------------------------------------
function Qvalue=Quant_Value(err);
% *************************************************************************
% This function is used as linear quantizer.The quantizer has totally
% 13 quantization level :
% 0,-+7,-+16,-+27,-+38,-+51,-+66
% *************************************************************************

if (abs(err)<=3) lev=0;
else if(abs(err)<=11) lev=7;
else if(abs(err)<=21) lev=16;
else if(abs(err)<=32) lev=27;
else if(abs(err)<=44) lev=38;
else if(abs(err)<=58) lev=51;
else lev=66;
end;
end;
end;
end;
end;
end;
Qvalue = lev;
if (err ~= 0)
Qvalue = lev*(err / abs(err));
end
%-------------------------------------------------------------------------
function Rvalue=Restor_Value(quan_err,pre_val)
% *************************************************************************
% This function is used to get restored value of DPCM
% x=^x+Quant_Error
% *************************************************************************

Rvalue=quan_err + pre_val;

%-------------------------------------------------------------------------
function Cvalue=Clip_Value(res_val)
% *************************************************************************
% This function is used to clip to restored value to 8_bit value
% 0 if x<0
% x'=255 if x>255
% x otherwise
% *************************************************************************

if (res_val<0) Cvalue=0;
else if(res_val>255) Cvalue=255;
else Cvalue=res_val;
end;
end;

四、 实验结果及分析（原图像和处理后的图像比较及分析）

如上图所示，即为实验所得的结果图像。由实验结果可以看出，一维和二维预测编码图像非常接近原图像；从误码图像中可以看出，一维预测编码会将误差延续到行末端，而二维预测编码其将误码延续至其后的斜后方，所以二维预测编码的误差会偏小。

实验名称：JPEG压缩编码

九、 实验目的
（1） 掌握nxn子块的DCT图像变换及频谱特点。
（2） 熟悉JPEG基本系统的图像编解码方法。

二、 实验内容
（1） 编程实现nxn子块DCT变换的图像频谱显示，8x8子块DCT变换系数按“Z”扫描图像重建，计算图像的均方根误差RMSE，显示误差图像和误差直方图。
（2） 编程实现JPEG压缩编码，进行8x8子块的DCT图像变换，JPEG量化矩阵的量化与反量化，8x8子块DCT的图像重建，计算图像的均方根误差RMSE，显示误差图像和误差直方图。

三、 实验方法及编程
用Matlab进行编程，主要是对照着实验书上面的算法流程图进行代码的编 写，主要的代码如下所示：
function newbuf=DctBlock(oldbuf,Block)
% *************************************************************************
% 函数名称：
% DctBlock() DCT n*n 块频谱函数
%
% 参数：
% oldbuf 原图像数组
% Block DCT n*n 当前块选择值
% newbuf 存放处理后的图像二维数组
%
% 说明：
% 根据Block块的当前选择值，计算原图像的n*n块DCT变换，并转换为可视频谱图，
% 有利于频谱的观察。
% *************************************************************************

oldbuf=double(oldbuf);% 添加程序
H=dctmtx(Block);
newbuf=blkproc(oldbuf,[Block Block],'P1*x*P2',H,H');
newbuf=log(abs(newbuf));
subplot(2,2,2);
imshow(newbuf,[]);

%-------------------------------------------------------------------------
function newbuf=DctCode(oldbuf,DCTch)
% *************************************************************************
% 函数名称：
% DctCode() DCT 8*8 块系数“Z”字扫描图像压缩函数
%
% 参数：
% oldbuf 原图像数组
% DCTch DCT 8*8 块“Z”扫描当前系数选择值
% newbuf 存放处理后的图像二维数组
%
% 说明：
% 计算图像的8×8子块DCT变换，按“Z”字扫描顺序，根据DCTch参数，只保留64个
% DCT系数中的前DCTch个系数，对修改后的DCT系数用逆DCT变换重建图像，得到DCT变
% 换的压缩图像。计算重建图像的均方根误差RMSE ；显示误差图像和误差直方图。
% *************************************************************************

zigzag=[1 2 6 7 15 16 28 29 %设置z扫描顺序
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];

tbuf=ones(8); %定义8*8全1数组

maskbuf=tbuf .* zigzag<=DCTch; %根据当前DCTch值得到“Z”字扫描的二值掩模
oldbuf=double(oldbuf);
H=dctmtx(8); % 添加程序
dctno=blkproc(oldbuf,[8 8],'P1*x*P2',H,H');
newbuf=blkproc(dctno,[8 8],'P1*(x.*P2)*P3',H',maskbuf,H);
subplot(2,2,2);
imshow(newbuf,[]);
eimag(oldbuf,newbuf);

%-------------------------------------------------------------------------
function newbuf=JpegCode(oldbuf)
% *************************************************************************
% 函数名称：
% JpegCode() JPEG近似基准编码
%
% 参数：
% oldbuf 原图像数组
% newbuf 存放处理后的图像二维数组
%
% 说明：
% 实现JPEG压缩编码，进行 8×8 子块DCT变换、JPEG量化矩阵的量化与反量化，
% 8×8子块DCT的图像重建；计算重建图像的均方根误差RMSE ；显示误差图像和误差直方图。
% *************************************************************************
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] ; %标准亮度量化表

oldbuf=double(oldbuf); %原图像数据转换为双精度
H=dctmtx(8); % 添加程序
dctno=blkproc(oldbuf,[8 8],'P1*x*P2',H,H');
jpegno=blkproc(dctno,[8 8],'round(x./P1)',z);
jpegno1=blkproc(jpegno,[8 8],'x.*P1',z);
newbuf=blkproc(jpegno1,[8 8],'P1*x*P2',H',H);
subplot(2,2,2);
imshow(newbuf,[]);
eimag(oldbuf,newbuf);

%-------------------------------------------------------------------------
function eimag(oldbuf,newbuf)
% *************************************************************************
% 函数名称：
% eimag()
% 说明：
% 计算重建图像的均方根误差RMSE ；显示误差图像和误差直方图。
% *************************************************************************
e = double(oldbuf)-newbuf; %计算原图像与压缩图像之差
[m,n]=size(e);RMSE=sqrt(sum(e(:).^2) / (m*n)); %计算均方根误差(root-mean-square error)
if RMSE %如果有误差，即rmse不为0
emax = max(abs(e(:))); %找图像差最大值
[h,x] = hist(e(:),emax); %用于生成直方图数据
if length(h) >= 1 %如果有数据
s=max(h(:));
subplot(2,2,3);
%imshow(hist);
bar(x,h/s,'k'); %显示图像差值直方图
RMSE = num2str(RMSE); %把数值转换为字符串
strRMSE = strcat('图像差值直方图 均方根误差RMSE= ',RMSE); %把多个串连接成长串
title(strRMSE);
e = mat2gray(e,[-emax,emax]); %显示差值图像
subplot(2,2,4),imshow(e);
title('原图像与压缩图像的差值图像'); %显示结果图像标题
end
end

四、 实验结果及分析（原图像和处理后的图像比较及分析）

如上图所示，即为实验的结果图像。对于DCT变换，变换后的高频部分存在于图像的左上角，代表着图像的轮廓，而低频部分代表图像的细节部分。对DCT系数做反DCT变换则可复原出原图。若反变换前对DCT系数进行取舍则可以降低码率，但是，会对图像质量带来一定的影响。如上图所示，当保留的系数越少，恢复的图像方块效应越严重，误差越大，差值图像越明显；当保留的系数越多，图像越接近于原图，误差越小，差值图像越不明显。