a simple dpcmcode.m

来自「本人自己编写的一个图像预测编码（DPCM）的简单程序。」· M 代码 · 共 76 行

%A DPCM COMPRESSION CODE
%Using format:X(n)=a1*X(n-1)+a2*X(n-2)
clc
load woman;
[m,n]=size(X);
r1=0;r2=0;r3=0;
for i=1:m
    r1=r1+X(i,1)^2+X(i,2)^2;
    r2=r2+X(i,1)*X(i,2);
    for j=3:n
        r1=r1+X(i,j)^2;
        r2=r2+X(i,j)*X(i,j-1);
        r3=r3+X(i,j)*X(i,j-2);
    end
end
r1=r1/(m*n);r2=r2/(m*n);r3=r3/(m*n);
a1=r2*(r1-r3)/(r1^2-r2^2);
a2=(r3*r1-r2^2)/(r1^2-r2^2);
Y=zeros(m,n);
c=zeros(m,n);
for i=1:m
    for j=3:n
        Y(i,j)=a1*X(i,j-1)+a2*X(i,j-2);
        c(i,j)=X(i,j)-Y(i,j);
    end
end
for i=2:m
    Y(i,1)=a1*X(i-1,1)+a2*X(i-1,2);
    c(i,1)=X(i,1)-Y(i,1);
    Y(i,2)=a1*X(i-1,2)+a2*X(i-1,1);
    c(i,2)=X(i,2)-Y(i,2);
end
c(1,1)=0;c(1,2)=0;
rc=abs(c);
rc=reshape(rc,m*n,1);
rc=sort(rc);
bit=8;
for i=1:(2^(bit-1))
    q=round(m*n*i/(2^(bit-1)));
    d(i)=rc(q);
end
e=zeros(m,n);
for i=2:(2^(bit-1))
    [k1,k2]=find(abs(c)<=d(i)&abs(c)>d(i-1));
    q1=length(k1);
    for j=1:q1
        if(c(k1(j),k2(j))>=0) e(k1(j),k2(j))=(d(i)+d(i-1))/2;
        else e(k1(j),k2(j))=-(d(i)+d(i-1))/2;
        end
    end
end
[k1,k2]=find(abs(c)<=d(1));
q=length(k1);
for j=1:q
    if(c(k1(j),k2(j))>=0) e(k1(j),k2(j))=d(1)/2;
    else e(k1(j),k2(j))=-d(1)/2;
    end
end
e(1,1)=X(1,1);
e(1,2)=X(1,2);
e=round(e*(2^(bit-1)))/(2^(bit-1));
%decode
cX=zeros(m,n);
cX(1,1)=e(1,1);
cX(1,2)=e(1,2);
for i=1:m
    if(i~=1)
      cX(i,1)=e(i,1)+a1*X(i-1,1)+a2*X(i-1,2);
      cX(i,2)=e(i,2)+a1*X(i-1,2)+a2*X(i-1,1);
    end
    for j=3:n
        cX(i,j)=e(i,j)+a1*cX(i,j-1)+a2*cX(i,j-2);
    end
end
subplot(221);image(X);colormap(map);title('原始图像');
subplot(222);image(cX);colormap(map);title('压缩图像');