rand('seed', 0)

im = im2double(rgb2gray(imread('lenna.png')));
[h w] = size(im);
len = numel(im);

x1 = im(:, 1:2:end);
x2 = im(:, 2:2:end);
X = [x1(:) x2(:)]';



% [U, J, I] = unique(X', 'rows');
K = 100;
[idx, C] = kmeans(X', K, 'replicates', 1, 'emptyaction', 'drop');
% idx = idx(I);
C = double(im2uint8(C));
% To transmit:
% Centres, C
% Codes, idx (nearest centres)

% Reconstruct the image from codes
Xq = C(idx, :)';
imrec = zeros(size(im));
imrec(:, 1:2:end) = reshape(Xq(1, :), h, w/2);
imrec(:, 2:2:end) = reshape(Xq(2, :), h, w/2);
imrec = imrec ./ 255;

figure(4), imshow(im)
figure(5), imshow(imrec)
figure(6), imsc(abs(im - imrec))

lenh = len_huffman(im(:));
lena = len_arith(im(:));
fprintf('Naively:\n');
fprintf('Original: %d,\n Huffman: %d (Ratio: %f),\n Arithmetic: %d (Ratio: %f)\n', ...
    len*8, lenh, lenh/(len*8), lena, lena/(len*8));

trans = [C(:); idx(:)];
lenhq = len_huffman(trans);
lenaq = len_arith(trans);
fprintf('After quantisation:\n');
fprintf('Original: %d,\n Huffman: %d (Ratio: %f),\n Arithmetic: %d (Ratio: %f)\n', ...
    len*8, lenhq, lenhq/(len*8), lenaq, lenaq/(len*8));

