%用 J(k,l)=I(k,l)-(I(k,l-1)/2+I(k-1,l)/4+I(k-1,l-1)/8+I(k-1,l+1)/8)進行預測編碼 
clc,clear all
I = rgb2gray(imread('lena.jpg')); %彩色圖像轉(zhuǎn)換為灰度圖像
I2=I; 
I=double(I); 
fid = fopen('mydata.dat','w'); 
[m n] = size(I); 
J=zeros(m,n); 
%邊界點不可預測，保留原值
J(1:m,1)=I(1:m,1); 
J(1,1:n)=I(1,1:n); 
J(1:m,n)=J(1:m,n); 
J(m,1:n)=J(m,1:n);
%計算預測器系數(shù)
J1=zeros(1,m*n);
R1=zeros(1,5);
A1=zeros(1,4); %預測器系數(shù)
B1=zeros(4,4);
B2=zeros(4,1);
for k =1:m
    for l=1:n
        J1((k-1)*m+l)=I(k,l);
    end
end
for k=0:4
    for l=k+1:20
        R1(k+1)=R1(k+1)+J1(1,l-k)*J1(1,l)/20;
    end
end
for k=1:4
    for l=1:4
        B1(k,l)=R1(abs(k-l)+1);
    end
end
for k=1:4
    B2(k,1)=R1(k+1);
end
A1=B1\B2;
A1
%計算預測誤差
for k=2:m-1 
    for l=2:n-1 
        J(k,l)=I(k,l)-(I(k,l-1)*A1(1)+I(k-1,l)*A1(2)+I(k-1,l-1)*A1(3)+I(k-1,l+1)*A1(4)); 
    end 
end 
%量化預測誤差
a=max(J);
b=min(J);
c=20; %量化層數(shù)  
t=round((a+b)/c);  %量化間隔
for k=a:t:b
    for p=1:m
        for q=1:n
            if((J(p,q)>k) & (J(p,q)<=(k+t)))
                J(p,q)=(2*k+t)/2;
            end
        end
    end
end
J=round(J);
%保存預測編碼文件
cout=fwrite(fid,J,'int8'); 
cc=fclose(fid); 
%預測編碼的解壓縮過程
fid=fopen('mydata.dat','r');
I1=fread(fid,inf,'int8');
tt=1;
for l=1:n
    for k=1:m
        I(k,l) = I1(tt);
        tt=tt+1;
    end
end
I=double(I);
J=zeros(m,n);
J(1:m,1)=I(1:m,1);
J(1,1:n)=I(1,1:n);
J(1:m,n)=I(1:m,n);
J(m,1:n)=I(m,1:n);
%由預測誤差還原圖像
for k=2:m-1
    for l=2:n-1
        J(k,l)=I(k,l)+(J(k,l-1)*A1(1)+J(k-1,l)*A1(2)+J(k-1,l-1)*A1(3)+J(k-1,l+1)*A1(4));
    end
end
cc=fclose(fid);
J=uint8(J);
sqnr=20*log10(norm(double(I2))/norm(double(I2)-double(J)))
subplot(1,2,1),imshow(I2);
subplot(1,2,2),imshow(J);