close all; clear all;

image = imread('neptune.jpg');

R = image(:,:,1);
G = image(:,:,2);
B = image(:,:,3);

yCbCrImage = rgb2ycbcr(image);

Y = yCbCrImage(:,:,1);
Cb = yCbCrImage(:,:,2);
Cr = yCbCrImage(:,:,3);

sizeY = size(Y);
numCols = sizeY(2);
numRows = sizeY(1);

% Downsampe chrominance - by 4, horizontally and vertically
Cb_4_1_1 = downsample(Cb',4)';
Cr_4_1_1 = downsample(Cr',4)';
Cb_4_1_0 = downsample(Cb_4_1_1,4);
Cr_4_1_0 = downsample(Cr_4_1_1,4);

Cb_4_1_0 = imresize(Cb_4_1_0,[numRows numCols],'nearest');
Cr_4_1_0 = imresize(Cr_4_1_0,[numRows numCols],'nearest');

ximage_4_1_0(:,:,1) = Y;
ximage_4_1_0(:,:,2) = Cb_4_1_0;
ximage_4_1_0(:,:,3) = Cr_4_1_0;

%DCT
dct_ximage_4_1_0(:,:,1) = dct2(ximage_4_1_0(:,:,1));
dct_ximage_4_1_0(:,:,2) = dct2(ximage_4_1_0(:,:,2));
dct_ximage_4_1_0(:,:,3) = dct2(ximage_4_1_0(:,:,3));

%Remove upper 3/4
smallDCT = dct_ximage_4_1_0(1:numRows/4,1:numCols/4,1:3);
smallDCT_out = uint8(smallDCT);

%Quantize
trainingData(1:1024,1:16,1:3) = zeros(1:1024,1:16,1:3);
a = 1;
for i = 1:4:128
    for j = 1:4:128
        temp = reshape(dct_ximage_4_1_0(i:(i+3),j:(j+3),1:3),1,16,3);
        trainingData(a,:,:) = temp;
        a=a+1;
    end
end

[IDX1,C1,SUMD1,D1] = kmeans(trainingData(:,:,1),16);
[IDX2,C2,SUMD1,D1] = kmeans(trainingData(:,:,2),16);
[IDX3,C3,SUMD1,D1] = kmeans(trainingData(:,:,3),16);

reconstruct(1:1024,1:16,1:3) = zeros(1:1024,1:16,1:3);
for i = 1:1024
    reconstruct(i,:,1) = C1(IDX1(i,:),:);
    reconstruct(i,:,2) = C2(IDX2(i,:),:);
    reconstruct(i,:,3) = C3(IDX3(i,:),:);
end

a = 1;
reconstructed_ximage_4_1_0 = zeros(1:128,1:128,1:3);
for i = 1:4:128
    for j = 1:4:128
        temp = reshape(reconstruct(a,1:16,1:3),4,4,3);
        reconstructed_ximage_4_4_2(i:(i+3),j:(j+3),1:3) = temp;
        a = a + 1;
    end
end

%THIS CONSTITUTES COMPRESSED IMAGE
sI1 = rle(IDX1);
sI2 = rle(IDX2);
sI3 = rle(IDX3);

[nc sizeIdx1] = size(sI1{1});
[nc sizeIdx2] = size(sI2{1});
[nc sizeIdx3] = size(sI3{1});

totalBits = (sizeIdx1 + sizeIdx2 + sizeIdx3)*4 + 16*16*3*32;

pic(:,:,1) = idct2(smallDCT(:,:,1),512,512);
pic(:,:,2) = idct2(smallDCT(:,:,2),512,512);
pic(:,:,3) = idct2(smallDCT(:,:,3),512,512);

imwrite(ycbcr2rgb(uint8(pic)),'Restored_Neptune_a.jpg');
