function xC = calic(A)
% calic           CALIC-like lossless coding of an image (integer matrix)
% This is a simplified variant of CALIC algorithm, actually only the
% (linerar) prediction and the context selection.
% It is assumed that the dynamic range of the image is 255. The image, i.e.
% prediction error, is mapped into sequences, depending on energy estimates. 
%
% The inverse is not ready yet, but is planned to be in this file as 
% Ar = calic(xC, size(A)), decode xC so that Ar == A.
% 
% ex:  A = double(imread('lena.bmp'))-128;
%      xC = calic(A);  [b,r] = estimateBits(xC, 'Huff06', 2);

% øving 7 i MIK230 2009 har mer om tapsfri koding av bilde
%    test comparing 'calic' to DPCM and mypred
% A = double(imread('lena.bmp')) - 128;
% r1 = imageapprox(A, 'd','ex312Jan191930.mat', 'tPSNR',38, 'v',1, 'del',-10);
%    -- kun DC komponent  (range not -128 to 128)
% dc = reshape(r1.Zdc,64,64);
% xCdc = calic(dc); 
% [b,r] = estimateBits(xCdc, 'Huff06', 2);   % 16904 bits
% xCdc = {[dc(1); diff(dc(:))]};             
% [b,r] = estimateBits(xCdc, 'Huff06', 2);   % 17520 bits
% xCdc = mypred(dc);
% [b,r] = estimateBits(xCdc, 'Huff06', 2);   % 16200 bits        !!
%    -- hele bildet
% xC = calic(A); 
% [b,r] = estimateBits(xC, 'Huff06', 2);   % bit rate 3.9550   !!!
% xC = mypred(A, 'nofS',11); 
% [b,r] = estimateBits(xC, 'Huff06', 2);   % bit rate 4.1594

visMye = 1;
nofS = 0;   % nofS = 0 is more simple and usually good
            %   nofS = 3; appropriate for dc above

if visMye
    disp('CALIC prediksjonsmetode.');
end
% Måten å dele opp i sekvenser for koding, Sekv, tilsvarer (grovt) context
% i CALIC. Jeg har tilpasset disse noe, de lages ikke i CALIC men er her
% hensiktsmessige siden entropikoding gjøres utenfor denne funksjonen, men
% er integrert i CALIC.

X = double(A);
[ROWS,COLS] = size(X);   % OBS NXM
Y = zeros(ROWS,COLS);
Sekv = ones(size(X));              % utgangspunkt sekv 1
if (nofS > 0); EE = 100*Sekv; end;
Y(1,1) = X(1,1);                   % første pixel ingen prediksjon
Y(2:ROWS,1) = X(2:ROWS,1)-X(1:(ROWS-1),1);  % venstre kolonne (prediksjon er pixel over)
j = 1;
for i=4:ROWS
    ee = 3*abs(X(i-3,j)-X(i-2,j))+5*abs(X(i-2,j)-X(i-1,j));
    if (nofS > 0); EE(i,j) = ee; end;
    Sekv(i,j) = ee2s(ee);
end
Y(1,2:COLS) = X(1,2:COLS)-X(1,1:(COLS-1));  % øverste linje (prediksjon er pixel til venstre)
i = 1;
for j=4:COLS
    ee = 3*abs(X(i,j-3)-X(i,j-2))+5*abs(X(i,j-2)-X(i,j-1));
    if (nofS > 0); EE(i,j) = ee; end;
    Sekv(i,j) = ee2s(ee);
end
%
for i=2:ROWS
    for j=2:COLS
        xW = X(i,j-1);        % pixel til venstre
        if j==2
            xWW = X(i,j-1);
        else
            xWW = X(i,j-2);   % pixel to hakk til venstre
        end
        xN = X(i-1,j);       % pixel over
        xNW = X(i-1,j-1);    % pixel over til venstre
        if i==2
            xNN = xN;
        else
            xNN = X(i-2,j);      % pixel to hakk over
        end
        if j==COLS
            xNE = xN;
            xNNE = xNN;
        else
            xNE = X(i-1,j+1);    % pixel over til høyre
            if i==2
                xNNE = xNE;
            else
                xNNE = X(i-2,j+1);   % pixel to hakk over og til høyre
            end
        end
        dh = abs(xW-xWW)+abs(xN-xNW)+abs(xNE-xN);
        dv = abs(xW-xNW)+abs(xN-xNN)+abs(xNE-xNNE);
        D = dv-dh;
        ee = dv+dh;  % energi estimat
        if (D > 80) % skarp horisontal kant
            p = xW;        % pixel til venstre
        elseif (D < -80) % skarp vertical kant
            p = xN;     % pixel over,  sekvens 1
        else
            p = (xN+xW)/2+(xNE-xNW)/4;   % prediction 2
            if (D > 32) % horisontal kant
                p = (p+xW)/2;
            elseif (D > 8)   % svak horisontal kant
                p = (3*p+xW)/4;
            elseif (D < -32)    % vertikal kant
                p = (p+xN)/2;
            elseif (D < -8)   % svak vertikal kant
                p = (3*p+xN)/4;
            end
            ee = ee+2*abs(p-xW);  % energi estimat (Delta) >= 0
        end
        if (nofS > 0); EE(i,j) = ee; end;
        Sekv(i,j) = ee2s(ee);
        %
        p = floor(p+0.4);   % avrundes med svak tendens nedover
        Y(i,j) = X(i,j)-p;
    end
end
%
if (nofS > 0)  % sekvenser velges på en mer komplisert måte
    t = sort(EE(:));
    vLim = [ceil(t( floor((1:(nofS-1))*((ROWS*COLS)/nofS)) ))', inf];
    disp(vLim);
    Sekv = zeros(ROWS,COLS);
    for i = 1:ROWS
        for j = 1:COLS
            Sekv(i,j) = find(EE(i,j) <= vLim, 1);
        end
    end
end

% danner sekvensene ut fra Sekv matrisa
Imax = max(Sekv(:));
xC = cell(Imax,1);
for i=1:Imax
    xC{i} = Y(Sekv==i);
    if visMye && (numel(xC{i}) > 0)
        disp(['Sekvens ',int2str(i),' har ',int2str(numel(xC{i})),' forekomster.',...
            '  min = ',int2str(min(xC{i})),'  max = ',int2str(max(xC{i})),...
            '  mean = ',num2str(mean(xC{i})), '  std = ',num2str(std(xC{i}))]);
    end
end
% For alle prediksjonsmetoder  har en at rekonstruksjon er greitt
% og gjøres ikke her, Xr = X;

return;


function s = ee2s(ee)
if ee > 140
    s = 1;
elseif ee > 85
    s = 2;
elseif ee > 60
    s = 3;
elseif ee > 42
    s = 4;
elseif ee > 25
    s = 5;
elseif ee > 15
    s = 6;
elseif ee > 5
    s = 7;
else
    s = 8;
end
return;
