clear all; clearvars; close all; clc;

% Displays B&W Photo
figure('Name','Magilla Gorilla');
x=imread('baboonbw.png');   % Reads in photo
imagesc(x);                 % Displays photo in Matlab figure
colormap(gray);

% Part 3a (Extract and display an 8x8 pixel block)
figure('Name','Magilla Gorilla ===> (8x8 Pixel Block:  Raw)');
X=x(81:88,81:88)
imagesc(X);                 % Displays 8x8 pixel block in figure
colormap(gray);

% Part 3b (Apply 2-Dimensional Discrete Cosine Transform (2D-DCT)
%   2D-DCT Matrix C
n=8;
C = zeros(n);
for i=1:n
    for j=1:n
        C(i,j)=cos((i-1)*(2*j-1)*pi/(2*n));
    end
end
C=sqrt(2/n)*C;
C(1,:)=C(1,:)/sqrt(2);
%C(1,:)=1/sqrt(2);
C_temp=C

% Quantization (Compression:  less important frequencies discarded)
p=1
Q = p*8./hilb(8);
X_double = double(X);
X_shiftd128 = X_double - 128

% Apply the 2D-DCT
%Y = dct(dct(X_shiftd128')')
Y=C*X_shiftd128*C';
Y_quant = round(Y./Q)

% Recover the image:
%   1.  De-quantization (Decompression)
%   2.  Inverse DCT
%   3.  Offset 128
%   4.  convert to uint8
Y_dequant = Y_quant.*Q
X_dequant = C'*Y_dequant*C;
X_offset128 = X_dequant + 128;
X_final = uint8(X_offset128);
figure('Name','Magilla Gorilla ===> (8x8 Pixel Blk:  Compressed/Decompressed)');
imagesc(X_final);
colormap(gray);
axis square;
X =

  8×8 uint8 matrix

   113   124   176   203   183   141   167   133
   190   183   169   157   137   158   144   126
   177   175   145   108    95    61    55    70
   158   140   113    98   158   179   197   199
   184   144   135   148   151   162   175   155
   190   190   201   206   210   203   186   175
   191   210   208   181   155   157   151   125
   143   125   136   142   123   153   148   176


C_temp =

  Columns 1 through 7

    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.4619    0.1913   -0.1913   -0.4619   -0.4619   -0.1913    0.1913
    0.4157   -0.0975   -0.4904   -0.2778    0.2778    0.4904    0.0975
    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.1913   -0.4619    0.4619   -0.1913   -0.1913    0.4619   -0.4619
    0.0975   -0.2778    0.4157   -0.4904    0.4904   -0.4157    0.2778

  Column 8

    0.3536
   -0.4904
    0.4619
   -0.4157
    0.3536
   -0.2778
    0.1913
   -0.0975


p =

     1


X_shiftd128 =

   -15    -4    48    75    55    13    39     5
    62    55    41    29     9    30    16    -2
    49    47    17   -20   -33   -67   -73   -58
    30    12   -15   -30    30    51    69    71
    56    16     7    20    23    34    47    27
    62    62    73    78    82    75    58    47
    63    82    80    53    27    29    23    -3
    15    -3     8    14    -5    25    20    48


Y_quant =

    28     3     0     0     0     0     0     0
    -4     1    -1     0     0     0     0     0
     0     0    -1    -1     0     0     0     0
     4    -1    -2     0     0     0     0     0
    -1    -3     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
    -1     0     0     0     0     0     0     0
    -1     1     0     0     0     0     0     0


Y_dequant =

   224    48     0     0     0     0     0     0
   -64    24   -32     0     0     0     0     0
     0     0   -40   -48     0     0     0     0
   128   -40   -96     0     0     0     0     0
   -40  -144     0     0     0     0     0     0
     0     0     0     0     0     0     0     0
   -56     0     0     0     0     0     0     0
   -64    72     0     0     0     0     0     0