clear;

%Read file into a matrix
array_input_rgb = imread('blue.jpg'); 

%Reads file info
array_info  = imfinfo('blue.jpg'); 
width=array_info.Width;
height=array_info.Height;

% Colour conversion from RGB to YCbCr
array_input=rgb2ycbcr(array_input_rgb);

% Create matrix which its height and width are divisible by 8.
% If is not, fill up with zeros
W=ceil(width/8);
H=ceil(height/8);
input=zeros(H*8,W*8,3,'uint8');
for I=1:3
    input(1:height,1:width,I)=array_input(1:height,1:width,I);
end

% Divide input matrix into WxH 8x8 matrices
X=zeros(H,W,8,8,3,'uint8');
for I=1:3
    for J=1:H    
        for K=1:W        
            for row=1:8            
                for col=1:8                
                    X(J,K,row,col,I)=input((J-1)*8+row,(K-1)*8+col,I);            
                end
            end
        end
    end
end

% Define luminance quantization matrix
QY=[...
    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];

% Define chrominance quantization matrix
QC=[...
    17 18 24 47 99 99 99 99
    18 21 26 66 99 99 99 99
    24 26 56 99 99 99 99 99
    47 66 99 99 99 99 99 99
    99 99 99 99 99 99 99 99
    99 99 99 99 99 99 99 99
    99 99 99 99 99 99 99 99
    99 99 99 99 99 99 99 99];

% Define matrix transformation
CF=[...
  0.3536  0.3536  0.3536  0.3536  0.3536  0.3536  0.3536  0.3536                                 
  0.4904  0.4157  0.2778  0.0975 -0.0975 -0.2778 -0.4157 -0.4904       
  0.4619  0.1913 -0.1913 -0.4619 -0.4619 -0.1913  0.1913  0.4619                             
  0.4157 -0.0975 -0.4904 -0.2778  0.2778  0.4904  0.0975 -0.4157
  0.3536 -0.3536 -0.3536  0.3536  0.3536 -0.3536 -0.3536  0.3536                       
  0.2778 -0.4904  0.0975  0.4157 -0.4157 -0.0975  0.4904 -0.2778                
  0.1913 -0.4619  0.4619 -0.1913 -0.1913  0.4619 -0.4619  0.1913                       
  0.0975 -0.2778  0.4157 -0.4904  0.4904 -0.4157  0.2778 -0.0975];

%   a  a  a  a  a  a  a  a                                 
%   b  c  f  g -g -f -c -b       
%   d  e -e -d -d -e  e  d                             
%   c -g -b -f  f  b  g -c
%   a -a -a  a  a -a -a  a                       
%   f -b  g  c -c -g  b -f                
%   e -d  d -e -e  d -d  e                       
%   g -f  c -b  b -c  f -g
  
a = 0.3536;
b = 0.4904;
c = 0.4157;
d = 0.4619;
e = 0.1913;
f = 0.2778;
g = 0.0975;

trunk=zeros(H*8,W*8,3);
for I=1:3
    for J=1:H
        for K=1:W        
            % Create temporary matrix 
            temp=zeros(8,8,'uint8');        
            dct_stage1=zeros(8,8,'double');   
            dct_stage2=zeros(8,8,'double');  
            idct_stage1=zeros(8,8,'double');   
            idct_stage2=zeros(8,8,'double'); 
            
            % Processing 8X8 matrix
            temp(:,:)=X(J,K,:,:,I);         
            
            % Convert to double 
            temp=double(temp);        
            
            % Shift mean to zero
            temp=temp-128;             
            
            % Transpose
            temp=temp';   
            
            % 1-D DCT Transform
            for(count=1:8)
                dct_stage1(1,count)=(a*(temp(1,count) + temp(8,count))) + (a*(temp(2,count) + temp(7,count))) + (a*(temp(3,count) + temp(6,count))) + (a*(temp(4,count) + temp(5,count)));
                dct_stage1(2,count)=(b*(temp(1,count) - temp(8,count))) + (c*(temp(2,count) - temp(7,count))) + (f*(temp(3,count) - temp(6,count))) + (g*(temp(4,count) - temp(5,count)));
                dct_stage1(3,count)=(d*(temp(1,count) + temp(8,count))) + (e*(temp(2,count) + temp(7,count))) - (e*(temp(3,count) + temp(6,count))) - (d*(temp(4,count) + temp(5,count)));
                dct_stage1(4,count)=(c*(temp(1,count) - temp(8,count))) - (g*(temp(2,count) - temp(7,count))) - (b*(temp(3,count) - temp(6,count))) - (f*(temp(4,count) - temp(5,count)));
                dct_stage1(5,count)=(a*(temp(1,count) + temp(8,count))) - (a*(temp(2,count) + temp(7,count))) - (a*(temp(3,count) + temp(6,count))) + (a*(temp(4,count) + temp(5,count)));
                dct_stage1(6,count)=(f*(temp(1,count) - temp(8,count))) - (b*(temp(2,count) - temp(7,count))) + (g*(temp(3,count) - temp(6,count))) + (c*(temp(4,count) - temp(5,count)));
                dct_stage1(7,count)=(e*(temp(1,count) + temp(8,count))) - (d*(temp(2,count) + temp(7,count))) + (d*(temp(3,count) + temp(6,count))) - (e*(temp(4,count) + temp(5,count)));
                dct_stage1(8,count)=(g*(temp(1,count) - temp(8,count))) - (f*(temp(2,count) - temp(7,count))) + (c*(temp(3,count) - temp(6,count))) - (b*(temp(4,count) - temp(5,count)));
            end           
            
            % Transpose
            temp=dct_stage1';  
            
            % 1-D DCT Transform
            for(count=1:8)
                dct_stage2(1,count)=(a*(temp(1,count) + temp(8,count))) + (a*(temp(2,count) + temp(7,count))) + (a*(temp(3,count) + temp(6,count))) + (a*(temp(4,count) + temp(5,count)));
                dct_stage2(2,count)=(b*(temp(1,count) - temp(8,count))) + (c*(temp(2,count) - temp(7,count))) + (f*(temp(3,count) - temp(6,count))) + (g*(temp(4,count) - temp(5,count)));
                dct_stage2(3,count)=(d*(temp(1,count) + temp(8,count))) + (e*(temp(2,count) + temp(7,count))) - (e*(temp(3,count) + temp(6,count))) - (d*(temp(4,count) + temp(5,count)));
                dct_stage2(4,count)=(c*(temp(1,count) - temp(8,count))) - (g*(temp(2,count) - temp(7,count))) - (b*(temp(3,count) - temp(6,count))) - (f*(temp(4,count) - temp(5,count)));
                dct_stage2(5,count)=(a*(temp(1,count) + temp(8,count))) - (a*(temp(2,count) + temp(7,count))) - (a*(temp(3,count) + temp(6,count))) + (a*(temp(4,count) + temp(5,count)));
                dct_stage2(6,count)=(f*(temp(1,count) - temp(8,count))) - (b*(temp(2,count) - temp(7,count))) + (g*(temp(3,count) - temp(6,count))) + (c*(temp(4,count) - temp(5,count)));
                dct_stage2(7,count)=(e*(temp(1,count) + temp(8,count))) - (d*(temp(2,count) + temp(7,count))) + (d*(temp(3,count) + temp(6,count))) - (e*(temp(4,count) + temp(5,count)));
                dct_stage2(8,count)=(g*(temp(1,count) - temp(8,count))) - (f*(temp(2,count) - temp(7,count))) + (c*(temp(3,count) - temp(6,count))) - (b*(temp(4,count) - temp(5,count)));
            end  
            
            % Quantization
            switch(I)
                case 1
                    for(count=1:8)
                        temp(1,count)=dct_stage2(1,count)/QY(1,count);
                        temp(2,count)=dct_stage2(2,count)/QY(2,count);
                        temp(3,count)=dct_stage2(3,count)/QY(3,count);
                        temp(4,count)=dct_stage2(4,count)/QY(4,count);
                        temp(5,count)=dct_stage2(5,count)/QY(5,count);
                        temp(6,count)=dct_stage2(6,count)/QY(6,count);
                        temp(7,count)=dct_stage2(7,count)/QY(7,count);
                        temp(8,count)=dct_stage2(8,count)/QY(8,count);
                    end
                case 2
                    for(count=1:8)
                        temp(1,count)=dct_stage2(1,count)/QC(1,count);
                        temp(2,count)=dct_stage2(2,count)/QC(2,count);
                        temp(3,count)=dct_stage2(3,count)/QC(3,count);
                        temp(4,count)=dct_stage2(4,count)/QC(4,count);
                        temp(5,count)=dct_stage2(5,count)/QC(5,count);
                        temp(6,count)=dct_stage2(6,count)/QC(6,count);
                        temp(7,count)=dct_stage2(7,count)/QC(7,count);
                        temp(8,count)=dct_stage2(8,count)/QC(8,count);
                    end
                case 3
                    for(count=1:8)
                        temp(1,count)=dct_stage2(1,count)/QC(1,count);
                        temp(2,count)=dct_stage2(2,count)/QC(2,count);
                        temp(3,count)=dct_stage2(3,count)/QC(3,count);
                        temp(4,count)=dct_stage2(4,count)/QC(4,count);
                        temp(5,count)=dct_stage2(5,count)/QC(5,count);
                        temp(6,count)=dct_stage2(6,count)/QC(6,count);
                        temp(7,count)=dct_stage2(7,count)/QC(7,count);
                        temp(8,count)=dct_stage2(8,count)/QC(8,count);
                    end
            end
          
            % Round off the quantized matrix
            temp=round(temp);     
            
            % Dequantization  
            switch(I)
                case 1
                    for(count=1:8)
                        temp(1,count)=temp(1,count)*QY(1,count);
                        temp(2,count)=temp(2,count)*QY(2,count);
                        temp(3,count)=temp(3,count)*QY(3,count);
                        temp(4,count)=temp(4,count)*QY(4,count);
                        temp(5,count)=temp(5,count)*QY(5,count);
                        temp(6,count)=temp(6,count)*QY(6,count);
                        temp(7,count)=temp(7,count)*QY(7,count);
                        temp(8,count)=temp(8,count)*QY(8,count);
                    end
                case 2
                    for(count=1:8)
                        temp(1,count)=temp(1,count)*QC(1,count);
                        temp(2,count)=temp(2,count)*QC(2,count);
                        temp(3,count)=temp(3,count)*QC(3,count);
                        temp(4,count)=temp(4,count)*QC(4,count);
                        temp(5,count)=temp(5,count)*QC(5,count);
                        temp(6,count)=temp(6,count)*QC(6,count);
                        temp(7,count)=temp(7,count)*QC(7,count);
                        temp(8,count)=temp(8,count)*QC(8,count);
                    end
                case 3
                    for(count=1:8)
                        temp(1,count)=temp(1,count)*QC(1,count);
                        temp(2,count)=temp(2,count)*QC(2,count);
                        temp(3,count)=temp(3,count)*QC(3,count);
                        temp(4,count)=temp(4,count)*QC(4,count);
                        temp(5,count)=temp(5,count)*QC(5,count);
                        temp(6,count)=temp(6,count)*QC(6,count);
                        temp(7,count)=temp(7,count)*QC(7,count);
                        temp(8,count)=temp(8,count)*QC(8,count);
                    end
            end
            
            % Transpose
            temp=temp';
            
            % 1-D IDCT Transform
            for(count=1:8)
                idct_stage1(count,1)=(a*temp(count,1))+(b*temp(count,2))+(d*temp(count,3))+(c*temp(count,4))+(a*temp(count,5))+(f*temp(count,6))+(e*temp(count,7))+(g*temp(count,8));
                idct_stage1(count,2)=(a*temp(count,1))+(c*temp(count,2))+(e*temp(count,3))-(g*temp(count,4))-(a*temp(count,5))-(b*temp(count,6))-(d*temp(count,7))-(f*temp(count,8));
                idct_stage1(count,3)=(a*temp(count,1))+(f*temp(count,2))-(e*temp(count,3))-(b*temp(count,4))-(a*temp(count,5))+(g*temp(count,6))+(d*temp(count,7))+(c*temp(count,8));
                idct_stage1(count,4)=(a*temp(count,1))+(g*temp(count,2))-(d*temp(count,3))-(f*temp(count,4))+(a*temp(count,5))+(c*temp(count,6))-(e*temp(count,7))-(b*temp(count,8));
                idct_stage1(count,5)=(a*temp(count,1))-(g*temp(count,2))-(d*temp(count,3))+(f*temp(count,4))+(a*temp(count,5))-(c*temp(count,6))-(e*temp(count,7))+(b*temp(count,8));
                idct_stage1(count,6)=(a*temp(count,1))-(f*temp(count,2))-(e*temp(count,3))+(b*temp(count,4))-(a*temp(count,5))-(g*temp(count,6))+(d*temp(count,7))-(c*temp(count,8));
                idct_stage1(count,7)=(a*temp(count,1))-(c*temp(count,2))+(e*temp(count,3))+(g*temp(count,4))-(a*temp(count,5))+(b*temp(count,6))-(d*temp(count,7))+(f*temp(count,8));
                idct_stage1(count,8)=(a*temp(count,1))-(b*temp(count,2))+(d*temp(count,3))-(c*temp(count,4))+(a*temp(count,5))-(f*temp(count,6))+(e*temp(count,7))-(g*temp(count,8));
            end
            
            % Transpose
            temp=idct_stage1';                
            
            % 1-D IDCT Transform
            for(count=1:8)
                idct_stage2(count,1)=(a*temp(count,1))+(b*temp(count,2))+(d*temp(count,3))+(c*temp(count,4))+(a*temp(count,5))+(f*temp(count,6))+(e*temp(count,7))+(g*temp(count,8));
                idct_stage2(count,2)=(a*temp(count,1))+(c*temp(count,2))+(e*temp(count,3))-(g*temp(count,4))-(a*temp(count,5))-(b*temp(count,6))-(d*temp(count,7))-(f*temp(count,8));
                idct_stage2(count,3)=(a*temp(count,1))+(f*temp(count,2))-(e*temp(count,3))-(b*temp(count,4))-(a*temp(count,5))+(g*temp(count,6))+(d*temp(count,7))+(c*temp(count,8));
                idct_stage2(count,4)=(a*temp(count,1))+(g*temp(count,2))-(d*temp(count,3))-(f*temp(count,4))+(a*temp(count,5))+(c*temp(count,6))-(e*temp(count,7))-(b*temp(count,8));
                idct_stage2(count,5)=(a*temp(count,1))-(g*temp(count,2))-(d*temp(count,3))+(f*temp(count,4))+(a*temp(count,5))-(c*temp(count,6))-(e*temp(count,7))+(b*temp(count,8));
                idct_stage2(count,6)=(a*temp(count,1))-(f*temp(count,2))-(e*temp(count,3))+(b*temp(count,4))-(a*temp(count,5))-(g*temp(count,6))+(d*temp(count,7))-(c*temp(count,8));
                idct_stage2(count,7)=(a*temp(count,1))-(c*temp(count,2))+(e*temp(count,3))+(g*temp(count,4))-(a*temp(count,5))+(b*temp(count,6))-(d*temp(count,7))+(f*temp(count,8));
                idct_stage2(count,8)=(a*temp(count,1))-(b*temp(count,2))+(d*temp(count,3))-(c*temp(count,4))+(a*temp(count,5))-(f*temp(count,6))+(e*temp(count,7))-(g*temp(count,8));
            end
            
            % Round off the matrix
            temp=round(idct_stage2);        
            temp=temp+128;        
            
            for row=1:8            
                for col=1:8                
                    trunk((J-1)*8+row,(K-1)*8+col,I)=temp(row,col);            
                end
            end
            
        end
    end
end

array_output=zeros(height,width,3);
array_output(:,:,:)=trunk(1:height,1:width,:);
array_output=uint8(array_output);

% Colour conversion from YCbCr to RGB 
array_output_rgb=ycbcr2rgb(array_output);

subplot(2,2,1), imshow(uint8(array_input_rgb)), title('Original Image');
subplot(2,2,2), imshow(uint8(array_output_rgb)), title('Compressed Image');

% MSE Calculation
count = 0;
for I=1:3
    for J=1:height
        for K=1:width
            reg =(double(array_output_rgb(J,K,I))-double(array_input_rgb(J,K,I)))^2;
            count = count + reg;
        end
    end
end
MSE = count/(width*height*3)
PSNR = 10*log10((255)^2/MSE)

imwrite(array_output_rgb,'bluenew.jpg','jpg');
