% dct_demo.m

dat = imread('test_image.tif');
[m,n] = size(dat(:,:,1));
m = floor(m/8)*8; % to make things easier: make dimensions multiples of 8
n = floor(n/8)*8;

r = dat(1:m,1:n,1);
g = dat(1:m,1:n,2);
b = dat(1:m,1:n,3);
Y  = uint8(      0.299   *r +  0.587   *g  +  0.114   *b);

figure();
imagesc(Y); title('original'); colormap(gray);
imwrite(Y,'original_bw.tif','TIFF','Compression','none');

% quantization matrix
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];

Ytransformed = zeros(m,n);

% loop over 8x8 blocks
for bx = 1:m/8
for by = 1:n/8
   block = double(Y(bx*8-7:bx*8, by*8-7:by*8))-128; % get block, and shift it
   block = dct2(block);                             % discrete cosine transform
   block = round(block ./ Q);                       % quantization
   Ytransformed(bx*8-7:bx*8, by*8-7:by*8) = block;
end
end


clear r g b Y
Y = uint8(zeros(m,n));

% reconstruct Y channel
for bx = 1:m/8
for by = 1:n/8
   block = Ytransformed(bx*8-7:bx*8, by*8-7:by*8); % get block
   block = block.*Q;                               % dequantization
   block = idct2(block);                           % inverse discrete cosine transform
   block = uint8(block + 128);                     % shift
   Y(bx*8-7:bx*8, by*8-7:by*8) = block;
end
end

figure();
imagesc(Y); title('reconstructed'); colormap(gray);
imwrite(Y,'reconstructed_bw.tif','TIFF','Compression','none');