数字图像处理 实验五：差分脉冲编码调制(DPCM)

DIP实验5：DPCM

实验目的

掌握帧内DPCM的编解码方法（预测、量化）及其压缩性能。

实验内容

1）编写一维DPCM编解码程序；
2）编写二维DPCM编解码程序（预测系数[1/2,1/8,1/4,1/8]）；
3）分别计算上两种情况中重建图像的PSNR。其中，重建图像和原图像的差值用绝对误差表示，可设置放大因子n以便观察误差分布情况。量化器采用15个量化分层，量化电平分别为（0， 5， 10， 17， 28， 39， 52， 67）。实验图像为cla0。
4）当人为加入单个误码，可观察到怎样的误码图案？

输出图像排列格式如下：

原图像	一维DPCM重建图像	二维DPCM重建图像
	一维DPCM误码图案	二维DPCM误码图案

需要编写的函数：
1）发送端：线性预测编解码器函数
function [newbuf,Qerror]=Dpcm_code(oldbuf,M,N,dim)
2）求预测值函数
function Pvalue=Predict_Value(row,col,N,dim)
3）求量化值函数
function Qvalue=Quant_Value(err)
4）求重建值函数
function Rvalue=Restor_Value(quan_err,pre_val)
5）限幅函数
function Cvalue=Clip_Value(res_val)
6）接收端：传输误码解码器函数
function newbuf=Error_Code(M,N,dim)
7）峰值信噪比

附：计算预测值，预测公式如下：

1. 1_D DPCM:

if j=1（若像素在图像左边第1列，则预测值设为128）  	^x[i,j]=128           
if j>1                        					  	^x[i,j]=x'[i,j-1]     

2. 2_D DPCM:

if i=1,j=1			
	^x[i,j]=128           
if i=1,j>1			
	^x[i,j]=x'[ i,j-1]   
if i>1,j=1 or 255	
	^x[i,j]=x'[i-1,j ]   
if i>1,j>1	(Pirsch's predictor)	
	^x[i,j]=1/2x'[ i,j-1]+1/8x'[i-1,j-1]+1/4x'[i-1,j ]+1/8x'[i-1,j+1]

---

参考代码

% 预测编码
OI = imread('CLA1.bmp');
subplot(2,3,1);
imshow(OI);
title("原图像");

OI = double(OI);
[M,N] = size(OI);

[newImg1,Qerror1] = Dpcm_code(OI,M,N,1);
newImg1 = uint8(newImg1);
subplot(2,3,2);
imshow(newImg1);
title("一维DPCM重建图像")

% 计算峰值信噪比
psnr1 = calPsnr(OI,newImg1)

[newImg2,Qerror2] = Dpcm_code(OI,M,N,2);
newImg2 = uint8(newImg2);
subplot(2,3,3);
imshow(newImg2);
title("二维DPCM重建图像")

% 计算峰值信噪比
psnr2 = calPsnr(OI,newImg2)

newImg=Error_Code(Qerror1,M,N,1);
newImg = uint8(newImg);
subplot(2,3,5);
imshow(newImg);
title("一维DPCM误码图案")

newImg=Error_Code(Qerror2,M,N,2);
newImg = uint8(newImg);
subplot(2,3,6);
imshow(newImg);
title("二维DPCM误码图案")


%1）发送端：线性预测编解码器函数
function [newImg,Qerror]=Dpcm_code(oldbuf,M,N,dim)
    global newbuf;
    newbuf = zeros(M,N);  % 预测值矩阵，定义为全局变量，因为其他函数中需要用到
    Qerror = zeros(M,N);
    for i = 1:M
        for j = 1:N
            preVal = Predict_Value(i,j,N,dim); % 得到该像素点的预测值
            error = oldbuf(i,j) - preVal;   % 与原来的像素值相减得到预测误差
            Qerror(i,j) = Quant_Value(error); % 将预测误差量化
            rValue = Restor_Value(Qerror(i,j),preVal); % 用预测值加上量化后的误差
            newbuf(i,j) = Clip_Value(rValue); % 限幅函数
        end
    end
    newImg = newbuf;
end
            
            
% 2）求预测值函数
function Pvalue=Predict_Value(row,col,N,dim)
    global newbuf;
    % dim 为1时是一维预测，dim 为2时是二维预测
    switch dim
        case 1
            if col == 1 % 图像的第一列没有前一个预测像素值
                Pvalue = 128;
            else 
                Pvalue = newbuf(row,col-1);
            end
        case 2
            if row == 1 && col == 1 % （1，1）像素
                Pvalue = 128;
            end
            if row == 1 && col > 1 % 第一行，列数大于1的像素，它的预测值只有前一个像素 
                Pvalue = newbuf(row,col-1);
            end
            if (row > 1 && col == 1) || (row > 1 && col == N)  % 第一列或最后一列，且行数大于1
                Pvalue = newbuf(row-1,col);
            end
            if row > 1 && col > 1 && col < N  % 其余情况，四个像素值加权
                Pvalue = (1/2)*newbuf(row,col-1) + (1/8)*newbuf(row-1,col-1) + (1/4)*newbuf(row-1,col) + (1/8)*newbuf(row-1,col+1);
            end
    end
end


% 3）求量化值函数
% 量化电平分别为（0，5，10，17，28，39，52，67） 范围对应（3,11,21,32,44,58,66）
function Qvalue=Quant_Value(err)
    % 判断量化电平
    if abs(err) <= 3
        level = 0;
    elseif abs(err) <= 11
        level = 5;
    elseif abs(err) <= 21
        level = 10;
    elseif abs(err) <= 32
        level = 17;
    elseif abs(err) <= 44
        level = 28;
    elseif abs(err) <= 58
        level = 39;
    elseif abs(err) <= 66
        level = 52;
    else
        level = 67;
    end
    if err ~= 0
        Qvalue  = level * (err/abs(err)); % 量化值的正负
    else 
        Qvalue = 0;
    end
end
   
    
% 4）求重建值函数
function Rvalue=Restor_Value(quan_err,pre_val)
    Rvalue = quan_err + pre_val;
end


% 5）限幅函数
function Cvalue=Clip_Value(res_val)
    if res_val < 0
        Cvalue = 0;
    elseif res_val > 255
            Cvalue = 255;
    else
            Cvalue = res_val;
    end
end


% 6）接收端：传输误码解码器函数
function newImg=Error_Code(qerror,M,N,dim)
    global newbuf;
    newbuf = zeros(M,N);
    qerror(50,50) = 100; % 将某一量化误差更改
    for i = 1:M
        for j = 1:N
            preVal = Predict_Value(i,j,N,dim);
            rValue = Restor_Value(qerror(i,j),preVal);
            newbuf(i,j) = Clip_Value(rValue);
        end
    end
    newImg = newbuf;
end


% 7）峰值信噪比
function PSNR = calPsnr(img1,img2)
    img1 = double(img1);
    img2 = double(img2);
    [M,N] = size(img1);
    % 峰值信噪比公式
    MES = sum(sum((img1-img2).^2))/(M*N);
    PSNR = 20 * log10(255/sqrt(MES));
end

实验结果