% 11.2 Computer Problem 5
clear all; clearvars; close all; clc;
x=imread('Sedona_AZ.JPG');
figure(1);
axis on;
imagesc(x);
title('Sedona, AZ');

% Convert to Grayscale
x=double(x);
r=x(:,:,1);
g=x(:,:,2);
b=x(:,:,3);
xgray=0.2126*r+0.7152*g+0.0722*b;
xgray=uint8(xgray);
%{
figure(2);
axis on;
imagesc(xgray);
colormap(gray)
title('Grayscale');
%}

% Discrete Cosine Transform Matrix C
n=8;
C = zeros(n);
for i=2:n
    for j=1:n
        C(i,j)=sqrt(2/n)*cos((i-1)*(j-1/2)*pi/n);
    end
end
C(1,1:n)=ones(1,n)/sqrt(n);

% Adjust image to an 8*8 pixel block based on Grayscale size
size8 = zeros(size(xgray)/8);
[row,col] = size(size8);
for i=1:row
    for j=1:col
        rb=r(i*8-7:i*8, j*8-7:j*8); %red
        gb=g(i*8-7:i*8, j*8-7:j*8); %green
        bb=b(i*8-7:i*8, j*8-7:j*8); %blue
        % Numerical Conversion; Centering
        rd=double(rb);
        rc=rd-128;
        gd=double(gb);
        gc=gd-128;
        bd=double(bb);
        bc=bd-128;
        % Apply the 2D-DCT
        %ry=dct(dct(rc')'); gy=dct(dct(gc')'); by=dct(dct(bc')');
        ry=C*rc*C';
        gy=C*gc*C';
        by=C*bc*C';
        % Quantization (Compression:  less important frequencies discarded)
        p=1;
        %p=2;
        %p=4;
        %p=16;
        Q=p*8./hilb(8);
        ryq=round(ry./Q);
        gyq=round(gy./Q);
        byq=round(by./Q);
        % Recover the image:
        %   1.  De-quantization (Decompression)
        %   2.  Inverse DCT
        %   3.  Offset 128
        %   4.  convert to uint8
        rydq=ryq.*Q;
        rxdq=C'*rydq*C;
        rxe=rxdq+128;
        gydq=gyq.*Q;
        gxdq=C'*gydq*C;
        gxe=gxdq+128;
        bydq=byq.*Q;
        bxdq=C'*bydq*C;
        bxe=bxdq+128;
        % Merge 8*8 pixel blocks
        Xe(i*8-7:i*8, j*8-7:j*8, 1)=rxe;
        Xe(i*8-7:i*8, j*8-7:j*8, 2)=gxe;
        Xe(i*8-7:i*8, j*8-7:j*8, 3)=bxe;
    end
end
% Image output
figure(3);
axis on;
imagesc(uint8(Xe));
title(['p = ',num2str(p)]);