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(:)]';

Q = sqrt(100);

X = X .* 255;
Xq = round(X ./ Q);
% To transmit:
% Quantization constant, Q
% Codes, Xq (nearest centres)

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

figure(1), imshow(im)
figure(2), imshow(imrec)
figure(3), 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 = Xq(:);
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));

