clear;
clc;
close all;
I=imread( 'lena.jpg');%讀取圖
I=rgb2gray(I);
figure;
subplot(2,2,1),imshow(I) ;
title('原始圖象');

I=double(I);

% 量化表  

mask = [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];
%求出8*8的DCT變換核T
n=8;
for i=1:n
    for j=1:n
        if (i==1)
           T(i,j)=sqrt(1/n);
        else
            T(i,j)=sqrt(2/n)*cos((i-1)*(2*j-1)*pi/(2*n));
        end
    end
end
% 分塊處理
for m=0:31  
    for n=0:31
        for i=1:8
            for j=1:8
                Bn(i,j)=I(8*m+i,8*n+j);
            end
        end
b=Bn;

%每一塊做DCT變換
b=T*b*T';
%利用量化表進行量化處理
b=round(b./mask);
% 每一塊Z字掃描
%c(1:64)={b(1,1),b(1,2),b(2,1),b(3,1),b(2,2),b(1,3),b(1,4),b(2,3),b(3,2),b(4,1),b(5,1),b(4,2),b(3,3),b(2,4),b(1,5),
   %      b(1,6),b(2,5),b(3,4),b(4,3),b(5,2),b(6,1),b(7,1),b(6,2),b(5,3),b(4,4),b(3,5),b(2,6),b(1,7),b(1,8),b(2,7),
      %   b(3,6),b(4,5),b(5,4),b(6,3),b(7,2),b(8,1),b(8,2),b(7,3),b(6,4),b(5,5),b(4,6),b(3,7),b(2,8),b(3,8),b(4,7),
       %  b(5,6),b(6,5),b(7,4),b(8,3),b(8,4),b(7,5),b(6,6),b(5,7),b(4,8),b(5,8,),b(6,7),b(7,6),b(8,5),b(8,6),b(7,7),
       %  b(6,8),b(7,8),b(8,7),b(8,8)}
c(1)=b(1,1) ;  
g=1;h=2;e=1; 
while g>=1&g<=8&h>=1&h<=8 
      e=e+1 ;
      c(e)=b(g,h) ;
      g=g+1 ;
      h=h-1 ;
end 
g=g-1 ;
h=h+1 ;
while g>=1&g<=8&h>=1&h<=8 
    if g>h 
        g=g+1 ;
        while g>=1&g<=8&h>=1&h<=8
            e=e+1; 
            c(e)=b(g,h); 
            g=g-1 ;
            h=h+1 ;
        end 
    else break 
    end 
    g=g+1 ;
    h=h-1 ;
    if g<h 
        h=h+1 ;
        while g>=1&g<=8&h>=1&h<=8 
            e=e+1; 
            c(e)=b(g,h); 
            g=g+1 ;
            h=h-1 ;
        end 
    else break 
    end 
    g=g-1 ;
    h=h+1 ;
end 
g=g-2 ;
h=h+1 ;
while g>=1&g<=8&h>=1&h<=8 
    if g>h
        h=h+1 ;
        while g>=1&g<=8&h>=1&h<=8 
            e=e+1; 
           c(e)=b(g,h) ;
           g=g-1 ;
           h=h+1 ;
        end 
    else break 
    end 
    g=g+1 ;
    h=h-1 ;
    if g<h 
        g=g+1 ;
        while g>=1&g<=8&h>=1&h<=8 
            e=e+1; 
            c(e)=b(g,h); 
            g=g+1 ;
            h=h-1 ;
        end 
    else break 
    end 
    g=g-1 ;
    h=h+1 ;
end 

%反Z掃描

D=[c(1)  c(2)  c(6)  c(7)  c(15) c(16) c(28) c(29)
   c(3)  c(5)  c(8)  c(14) c(17) c(27) c(30) c(43)
   c(4)  c(9)  c(13) c(18) c(26) c(31) c(42) c(44)
   c(10) c(12) c(19) c(25) c(32) c(41) c(45) c(54)
   c(11) c(20) c(24) c(33) c(40) c(46) c(53) c(55)
   c(21) c(23) c(34) c(39) c(47) c(52) c(56) c(61)
   c(22) c(35) c(38) c(48) c(51) c(57) c(60) c(62)
   c(36) c(37) c(49) c(50) c(58) c(59) c(63) c(64)];
D1=T'*D *T;           %反DCT變換
iq=D.*mask;           %反量化處理
iq=T'*iq*T;            %反量化后再反DCT處理

%反分塊處理
 for s=1:8   
     for k=1:8
         An(8*m+s,8*n+k)=D1(s,k);      %未反量化的反分塊
         IQ(8*m+s,8*n+k)=iq(s,k);        %經(jīng)過反量化的反分塊
     end
 end
    end
end
%顯示量化后的圖形
subplot(2,2,2),imshow(An,[-5,25]) ; 
title('量化后的圖象');

%計算量化后圖形的均方誤差
MSE1=sum(sum((I-An).^2))/(256*256) 
Q=255;
%并計算峰值信噪比
PSNR1=10*log10(Q*Q/MSE1)
%顯示經(jīng)過反量化處理的圖形
subplot(2,2,3),imshow(IQ,[0,255])   ;         
title('恢復的圖象');
%計算反量化后圖形的均方誤差
MSE2=sum(sum((I-IQ).^2))/(256*256) 
Q=255;
%并計算峰值信噪比
PSNR2=10*log10(Q*Q/MSE2)