20150604_JPEG_compression

close all
clear all
x = [ 52    55    61    66    70    61    64    73 ;
 63    59    55    90    109    85    69    72 ;
 62    59    68    113    144    104    66    73 ;
 63    58    71    122    154    106    70    69 ;
 67    61    68    104    126    88    68    70 ;
 79    65    60    70    77    68    58    75 ;
 85    71    64    59    55    61    65    83 ;
 87    79    69    68    65    76    78    94
 ];

% each entry in x falls in range [0,255] , so subtracting the midpoint 128 to modify the range to [-128,127]
 g = x-128*ones(size(x));

 N=8;
G=zeros(N,N,N,N);
GG=zeros(N,N);
alpha_u=1;
alpha_v=1;

%2D DCT
for u = 0:N-1
    for v = 0:N-1
        for x = 0:N-1
            for y = 0:N-1

                G(u+1,v+1,x+1,y+1) = 0.25* ((u~=0)+(u==0)/sqrt(2))* ((v~=0)+(v==0)/sqrt(2)) * g(x+1,y+1) * cos(pi*u*(2*x+1/2)/N) * cos(pi*v*(2*y+1/2)/N) ;

            end
        end
        GG(u+1,v+1) = sum(sum(G(u+1,v+1,:,:),3),4)
    end
end

%Quantization matrix
%This allows one to greatly reduce the amount of information in the high frequency components.
%This is done by simply dividing each component in the frequency domain by a constant for that component, and then rounding to the nearest integer

Q=[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]

% rounds each element of X to the nearest integer
B=round(GG./Q)


%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% REVERSE %%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%
XX = B.*Q

%N=8;
X=zeros(N,N,N,N);
x=zeros(N,N);

for i = 0:N-1
    for j = 0:N-1
        for u = 0:N-1
            for v = 0:N-1

                X(i+1,j+1,u+1,v+1) = 0.25* ((u~=0)+(u==0)/sqrt(2))* ((v~=0)+(v==0)/sqrt(2)) * XX(u+1,v+1) * cos(pi*u*(2*i+1/2)/N) * cos(pi*v*(2*j+1/2)/N) ;

            end
        end
        x(i+1,j+1) = sum(sum(X(i+1,j+1,:,:),3),4);
    end
end

x

uc = zeros(N,N);
uc = x+128

er = f - uc
ab_er = abs(er)
point_er = 0;

for i = 0:N-1
    for j = 0:N-1
        point_er = point_er + ab_er(i+1,j+1)
    end
end
point_er/64