function varargout = EMP8_icm(s)
% EMP8_icm         Do Iterated Conditional Modes (ICM) on a label image 
%                  EMP8: Energy Minimization using Potts model and 8 neighbors
%
% This function may use:
%   EMP8_icm_mex(..)      to do the actual work
%   EMP8_cost_mex(..)     to find cost for a labeling
%
% Use of function:
% ----------------
%   L = EMP8_icm(s)
%   [L,c] = EMP8_icm(s)
% 
%     L   is the resulting new labeling after ICM
%     c   is a struct with cost results in each field
%     s   is a struct as returned by EMP8_argin, see help for that file
% 
% Example:
% -------------------------------------------------------
%   clear all; tim = 4;  sigma = 5;  
%   s = EMP8_argin(tim, sigma, tim, 'v',2);   % use Ltrue as well
%   [L,c] = EMP8_icm(s);
%   % and plot results in a figure
%   figure(1); clf; 
%   figTitle = ['Labeling of test image ',int2str(tim),' after using EMP8-icm.'];
%   LIS_plotcol(reshape(L,s.I,s.J), figTitle, c.errorRate, -1, c.Etot, c.Edata, c.E4+c.E8);

%----------------------------------------------------------------------
% Copyright (c) 2013.  Karl Skretting.  All rights reserved.
% University of Stavanger, Signal Processing Group
% Mail:  karl.skretting@uis.no   Homepage:  http://www.ux.uis.no/~karlsk/
% 
% HISTORY:  dd.mm.yyyy
% Ver. 1.0  12.12.2013  Made file (very similar to EMP8_erode)
%----------------------------------------------------------------------

mfile = 'EMP8_icm';

if ( ~exist('EMP8_icm_mex','file') || ...
     ~exist('EMP8_cost_mex','file') )
    c = [mfile,': one or more needed mex-files are not availble.'];
    disp(c);
    if (nargout >= 1); varargout{1} = []; end;
    if (nargout >= 2); varargout{2} = c; end;
    return;
end

%% do the calulations
if isfield(s,'useRand')
    useRand = s.useRand;
else
    useRand = true;
end
param = [s.M,s.I,s.J,-1, s.bF4,s.bF8,s.swUpq,0, s.lambda, s.dR, s.minsumRi, s.maxsumRi, ...
         s.numPix1, s.numPix2, s.shapeFactor, useRand];   % 16 parameters

if ((s.bF4 < 0) || (s.bF8 < 0))
    % this should normally not happen
    disp([mfile,': Values are not assigned to bF4 or bF8 yet, do that.']);
    disp([' bF4=',int2str(s.bF4),', bF8=',int2str(s.bF8),', lambda=',num2str(s.lambda), ...
         ', bF8percent=',int2str(s.bF8percent)  ]);
    if (s.bF4 < 0) 
        param(6) = s.bF8percent; 
    else   % (s.bF8 < 0)  % and (bF4 >= 0)
        param(6) = floor(s.bF4*(s.bF8percent/100)+0.5);
    end
    if (s.swUpq == 1)
        Ci = EMP8_cost_mex(s.Ri, s.Li, param);
    elseif (s.swUpq == 2)
        Ci = EMP8_cost_mex(s.Ri, s.Li, s.Upq, param);
    elseif (s.swUpq == 3)
        Ci = EMP8_cost_mex(s.Ri, s.Li, s.Uim, s.Utab, param);
    end
    if (s.bF4 < 0)
        s.bF4 = Ci(5);
        s.bF8 = Ci(6);
    end
    if (s.dR < 0)
        s.dR = Ci(14); 
        s.minsumRi = Ci(15); 
        s.maxsumRi = Ci(16);
    end
    param = [s.M,s.I,s.J,-1, s.bF4,s.bF8,s.swUpq,0, s.lambda, s.dR, s.minsumRi, s.maxsumRi, ...
             s.numPix1, s.numPix2, s.shapeFactor, -1];   % 16 parameters
    if (s.verbose >= 1)
        disp([mfile,': The values assigned to bF4 or bF8 are']);
        disp([' bF4=',int2str(s.bF4),', bF8=',int2str(s.bF8),', lambda=',num2str(s.lambda), ...
            ', bF8percent=',int2str(s.bF8percent)  ]);
    end
elseif (s.verbose >= 1)
    disp(' ');
    disp([mfile,': some values are bF4 = ',int2str(s.bF4),', bF8 = ',int2str(s.bF8),...
          ', lambda = ',num2str(s.lambda),', bF8percent = ',int2str(s.bF8percent), ...
          ', swUpq = ',int2str(s.swUpq)]);
end

% do icm in mex-file, EMP8_icm_mex
tic;
I = s.I; 
J = s.J;
L = reshape(s.Li,I,J);
mask = [false(I,1), (L(:,1:(J-1)) ~= L(:,2:J))] | ...
    [(L(:,1:(J-1)) ~= L(:,2:J)), false(I,1)] | ...
    [false(1,J); (L(1:(I-1),:) ~= L(2:I,:))] | ...
    [(L(1:(I-1),:) ~= L(2:I,:)); false(1,J)] ;
param(4) = nnz(mask);   % parameter P (or s.P)
if (s.swUpq == 1)
    [L, c] = EMP8_icm_mex(s.Ri, s.Li, int32(find(mask)-1), param);
elseif (s.swUpq == 2)
    [L, c] = EMP8_icm_mex(s.Ri, s.Li, int32(find(mask)-1), s.Upq, param);
elseif (s.swUpq == 3)
    [L, c] = EMP8_icm_mex(s.Ri, s.Li, int32(find(mask)-1), s.Uim, s.Utab, param);
end
icmTime = toc;

icmUpdated = c(17);   % updated
if (s.verbose >= 1)
    disp([' ICM relabeled ',int2str(icmUpdated),' pixels.']);
end

%% return results
if (nargout >= 1)
    varargout{1} = L;  
end
if (nargout >= 2)
    updated = icmUpdated;   %  = nnz(L ~= s.Li);
    if (numel(s.Lt) == (s.I*s.J))
        er = nnz(L ~= s.Lt)/(s.I*s.J);
    else
        er = -1;
    end
    icmCost = sprintf(s.formatCost, 'Licm', c([10,13,11,12,9]), 100*er);
    if (er < 0) % it is nicer to print nothing rather the negative error rate
        icmCost = [icmCost(1:72),'.'];
    end

    if (s.verbose >= 1)
        disp(icmCost);
    end
    if exist('c','var')   % and isnumeric(c)
        % den normale situasjonen
        c(17:20) = [updated; icmTime; icmUpdated; er];
        c = costTilStruct(c);
    else
        disp([mfile,': probably nothing was done.']);
        if isfield(s,'Ci')
            c = s.Ci;
        else
            if (s.swUpq == 1)
                c = EMP8_cost_mex(s.Ri, L, param);
            elseif (s.swUpq == 2)
                c = EMP8_cost_mex(s.Ri, L, s.Upq, param);
            elseif (s.swUpq == 3)
                c = EMP8_cost_mex(s.Ri, L, s.Uim, s.Utab, param);
            end
            c(17:20) = [updated; icmTime; icmUpdated; er];
            c = costTilStruct(c);
        end
    end
    c.icmCost = icmCost; 
    c.icmP = param(4);
    c.icmUseRand = useRand;

    varargout{2} = c;
end

return



function s = costTilStruct(c)
if (numel(c)<20)
    c = [c(:);-ones(20-numel(c),1)];
end
s = struct('Ctot',c(1),'Cdata',c(2),'C4',c(3),'C8',c(4),...
           'bF4',c(5),'bF8',c(6),'n4',c(7),'n8',c(8),...
           'Etot',c(9),'Edata',c(10),'E4',c(11),'E8',c(12),...
           'lambda',c(13),'dR',c(14),'minsumRi',c(15),'maxsumRi',c(16),...
           'updated',c(17),'time',c(18),'icmUpdated',c(19),'errorRate',c(20));
return;
       
