function varargout = LIS_cost(Ri, Li, varargin)
% LIS_cost         Return (and show) the cost (energy) for a labeled image
%
%  Results are bot returned in scaled form, each term in range 0 to 1, as
%     Etot = Edata + lambda*(E4 + E8); 
%  or as (large) integers for sum in Ri (R scaled and represented as int16)
%     Ctot = Cdata + C4 + C8;
%  where Cdata is sum of the cost in Ri for labeling Li
%  and C4 = n4*bF4 and C8 = n8*bF8
%  and n4 is number of horizontal and vertical edges between labeled pixels
%  and n8 is number of diagonal 'edges' between labeled pixels
%
%   Use of function:
%   ----------------
%   res = LIS_cost(R, L, 'option',value, ...)
%   res = LIS_cost(R, L, Ltrue, 'option',value, ...)
% 
%     res  is a struct with results in each field
%     R is data cost matrix, size IxJxM, class of R should be (double)
%       or (int16) then it is interpreted as Ri (int16) of size MIJx1 (or MxIxJ)
%       or an integer (1-13) for a stored matrix, see LIS_getR.m
%     L is input label matrix, size IxJ, or Li (uint8) of size IJx1 (or IxJ)
%       0 may be used for minimizing data term, or m in 1:M for m*ones(I,J)
%     Ltrue is true labeling (if available), may be an integer 1-13, 
%       or it may be skipped as it is only used when displaying results
%     Additional arguments may be given as pairs: 'option',value, ...
%       (or options may be given in one (or more) struct or cell-array)
% 
%     Some options are:
%   'lambda'  the parameter lambda to use for Etot (and Ctot), default 2
%   'bF8percent'  how many percent of bF4 the bF8 value should be, default 70
%             this is used if bF4 is not given (bF4 < 0)
%   'bF4'     border factor for 4-neighbors pixels, default -1
%             if not given (< 0) the bF4 value is found from data and lambda
%   'bF8'     if bF4 is given bF8 should also be given
%   'FigNo'   number for figure window to display results in, see LIS_plotcol.m
%       default: 0 (no plot)
%   'FigText' text to use in figure (default 'Li');
%   'scaleRi' may be given, if not a suitable value is found if needed
%   'returnRi', 'returnLi', and 'returnLtrue' may be set to true to include 
%       these variables (long column vectors) in output struct
%   'dR', 'minsumRi' and 'maxsumRi' may also be given as input
%   'useMex' true or false (default) mex-implementation and m-file implementation 
%       should be equal to each other, and no significant speed difference. 
%   'v' or 'verbose' may be 0/false (default), 1/true or 2 (very verbose).
% 
%   Examples:
%   ---------
%   tim = 3;  
%   res1 = LIS_cost(tim, 0, 1, 'v',2 );
%   R = LIS_getR(tim);  [I,J,M] = size(R);  
%   [~, Lmin] = min(R,[],3); Ltrue = teximage(tim+200);
%   res2 = LIS_cost(R, Lmin, Ltrue, 'verbose',2 ); % res2 (almost) same as res1
%   sigma = 5; Lg = LIS_smooth(R, sigma);
%   res3 = LIS_cost(R, Lg, Ltrue, 'v',2, 'FigNo',1 ); 

mfile = 'LIS_cost';

%% Check if input arguments are given
if (nargin < 2)  % just check number of input arguments 
    t = [mfile,': arguments must be given, see help.'];
    disp(t);
    if nargout >= 1
        varargout{1} = struct('Error',t);
    end
    return
end

%% defaults, initial values
formatCost = 'Labeling %8s --> E = %7.5f + %5.2f*(%7.5f + %7.5f) = %7.5f, err. %6.3f%% ';
lambda = 2;
bF4 = -1;
bF8 = -1;
bF8percent = 70;
dR = -1;
minsumRi = -1;
maxsumRi = -1;
figno = 0;
verbose = 0;
returnRi = false;
returnLi = false;
returnLtrue = false;
useMex = false;
figText = 'Li';
    
%% get the options
if (nargin >= 3)
    nofOptions = nargin-2;  % options start at argument 3 
    optionNumber = 1;
    if isnumeric(varargin{optionNumber})  
        Ltrue = varargin{optionNumber};    % first option was Ltrue
        optionNumber = 2;       % next option
    end
end
fieldNumber = 1;
while (optionNumber <= nofOptions)
    if isstruct(varargin{optionNumber})
        sOptions = varargin{optionNumber}; 
        sNames = fieldnames(sOptions);
        opName = sNames{fieldNumber};
        opVal = sOptions.(opName);
        % next option is next field or next (pair of) arguments
        fieldNumber = fieldNumber + 1;  % next field
        if (fieldNumber > numel(sNames)) 
            fieldNumber = 1;
            optionNumber = optionNumber + 1;  % next pair of options
        end
    elseif iscell(varargin{optionNumber})
        sOptions = varargin{optionNumber}; 
        opName = sOptions{fieldNumber};
        opVal = sOptions{fieldNumber+1};
        % next option is next pair in cell or next (pair of) arguments
        fieldNumber = fieldNumber + 2;  % next pair in cell
        if (fieldNumber > numel(sOptions)) 
            fieldNumber = 1;
            optionNumber = optionNumber + 1;  % next pair of options
        end
    else
        opName = varargin{optionNumber};
        opVal = varargin{optionNumber+1};
        optionNumber = optionNumber + 2;  % next pair of options
    end
    % interpret opName and opVal
    if strcmpi(opName,'M')
        if isnumeric(opVal); M = floor(opVal(1)); end;
    end
    if strcmpi(opName,'I')
        if isnumeric(opVal); I = floor(opVal(1)); end;
    end
    if strcmpi(opName,'J')
        if isnumeric(opVal); J = floor(opVal(1)); end;
    end
    if strcmpi(opName,'scaleRi')
        if isnumeric(opVal); scaleRi = opVal(1); end;
    end
    if strcmpi(opName,'lambda')
        if isnumeric(opVal); lambda = opVal(1); end;
    end
    if strcmpi(opName,'bF4')
        if isnumeric(opVal); bF4 = opVal(1); end;
    end
    if strcmpi(opName,'bF8')
        if isnumeric(opVal); bF8 = opVal(1); end;
    end
    if strcmpi(opName,'bF8percent')
        if isnumeric(opVal); bF8percent = opVal(1); end;
    end
    if strcmpi(opName,'dR')
        if isnumeric(opVal); dR = opVal(1); end;
    end
    if strcmpi(opName,'minsumRi')
        if isnumeric(opVal); minsumRi = opVal(1); end;
    end
    if strcmpi(opName,'maxsumRi')
        if isnumeric(opVal); maxsumRi = opVal(1); end;
    end
    if strcmpi(opName,'FigNo')
        if isnumeric(opVal); figno = opVal(1); end;
    end
    if strcmpi(opName,'FigText')
        if ischar(opVal); figText = opVal; end;
    end
    if strcmpi(opName,'verbose') || strcmpi(opName,'v')
        if (islogical(opVal) && opVal); verbose = 1; end;
        if isnumeric(opVal); verbose = opVal(1); end;
    end
    if strcmpi(opName,'returnRi')
        if (islogical(opVal) && opVal); returnRi = true; end;
        if isnumeric(opVal); returnRi = (opVal(1) > 0); end;
    end
    if strcmpi(opName,'returnLi')
        if (islogical(opVal) && opVal); returnLi = true; end;
        if isnumeric(opVal); returnLi = (opVal(1) > 0); end;
    end
    if strcmpi(opName,'returnLtrue')
        if (islogical(opVal) && opVal); returnLtrue = true; end;
        if isnumeric(opVal); returnLtrue = (opVal(1) > 0); end;
    end
    if strcmpi(opName,'useMex') 
        if (islogical(opVal) && opVal); useMex = true; end;
        if isnumeric(opVal); useMex = (opVal(1) > 0); end;
    end
end

%% check input Ri, the data matrix
if (isnumeric(Ri) && (numel(Ri) == 1))
    Ri = LIS_getR(Ri);  
end
if ~( isa(Ri,'int16') || isa(Ri,'double') )
    disp([mfile,' first argument R is not int16 or double, this may cause problem.']);
end
if isa(Ri,'int16') 
    if size(Ri,2) > 1
        [M,I,J] = size(Ri);    % if not a column vector size should be MxIxJ here!
    end
    Ri = Ri(:);
else  % process as if double
    [I,J,M] = size(Ri);
    if ~exist('scaleRi','var')   % if this exist it was given as input parameter
        scaleRi = 32000/max(abs(Ri(:)));
    end
    Ri = int16(scaleRi*reshape(permute(Ri,[3,1,2]),M*I*J,1));   % and make it int16
end
% Ri should now be the (int16) column vector with M*I*J elements

%% check input L, the label matrix
if (isnumeric(Li) && (numel(Li) == 1))
    if (Li == 0)
%         if verLessThan('matlab','7.12')  % or somthing less ??
            [temp,Li] = min(reshape(Ri, M, I*J));   % integers in Ri may give a little
%         else
%             [~,Li] = min(reshape(Ri, M, I*J));   % integers in Ri may give a little
%         end
        % bit different labeling than using R (double) direcly from LIS_getR()
        Li = reshape(Li,I,J);
    else
        Li = Li*ones(I,J);
    end
end
if ~( isa(Li,'uint8') || isa(Li,'double') )
    disp([mfile,' second argument L is not uint8 or double, this may cause problem.']);
end
if exist('M','var')
    if ~((numel(Ri)/numel(Li)) == M)
        t = [mfile,': sizes of R and L do not match, see help.'];
        disp(t);
        if nargout >= 1
            varargout{1} = struct('Error',t);
        end
        return
    end
end
if size(Li,2) == 1
    if ~(exist('I','var') && exist('J','var'))
        t = [mfile,': label image size IxJ is not availble, see help.'];
        disp(t);
        if nargout >= 1
            varargout{1} = struct('Error',t);
        end
        return
    end
    Li = reshape(Li, I, J);
else
    [I,J] = size(Li);
end
M = floor(numel(Ri)/(I*J));
Li = uint8(Li(:));

%% check input Ltrue, the label matrix
if exist('Ltrue','var')
    if (isnumeric(Ltrue) && (numel(Ltrue) == 1))
        Ltrue = teximage(Ltrue+200);
    end
    Ltrue = uint8(Ltrue(:));
    if (numel(Ltrue) ~= (I*J))
        t = [mfile,': Ltrue is not expected size, ignore it.'];
        disp(t);
        clear Ltrue
    end
end

%% check mex-file
if ~exist('LIS_cost_mex','file')
    useMex = false;
end

%% Display info for input
if (verbose > 1)  % very verbose
    disp(' ');
    disp([mfile,' started ',datestr(now)]);
    disp([' Use Ri (',class(Ri),', ',int2str(size(Ri,1)),'x',int2str(size(Ri,2)),...
          'x',int2str(size(Ri,3)),') and Li (',class(Li),', ',int2str(size(Li,1)),...
          'x',int2str(size(Li,2)),') M=',int2str(M),', I=',int2str(I),', J=',int2str(J),'.']);
    if exist('Ltrue','var')
        disp([' We also have Ltrue (',class(Ltrue),', ',int2str(size(Ltrue,1)),...
            'x',int2str(size(Ltrue,2)),'), ',...
            'bF4=',int2str(bF4),', bF8=',int2str(bF8),', lambda=',num2str(lambda)]);
    end
end

%% do the calulations
param = [M,I,J,-1, bF4,bF8,-1,0, lambda,dR,minsumRi,maxsumRi];  % 12 parameters

if (bF4 < 0) 
    param(6) = bF8percent; 
elseif (bF8 < 0)  % and bF4 >= 0
    param(6) = floor(bF4*(bF8percent/100)+0.5);
end

if exist('Ltrue','var')
    if useMex
        Ctrue = LIS_cost_mex(Ri, Ltrue, param);
    else
        Ctrue = LIS_cost_local(Ri, Ltrue, param);
    end
    ttrue = sprintf(formatCost, 'Ltrue', Ctrue([10,13,11,12,9]), 0);
    if (verbose > 0)
        disp(ttrue);
    end
    if (bF4 < 0)
        param([5,6,10,11,12]) = Ctrue([5,6,14,15,16]);        
    end
end

if useMex
    Ci = LIS_cost_mex(Ri, Li, param);
else
    Ci = LIS_cost_local(Ri, Li, param);
end
if exist('Ltrue','var')
    er = nnz(Li ~= Ltrue)/(I*J);
else
    er = 0;
end
ti = sprintf(formatCost, 'Li', Ci([10,13,11,12,9]), 100*er);
if (verbose > 0)
    disp(ti);
end
if (figno > 0)
    figure(figno);
    cc = conncomp(reshape(Li,I,J));
    LIS_plotcol(reshape(Li,I,J), figText, er, ....
        cc.NumObjects, Ci(9), Ci(10), Ci(11)+Ci(12))
end


%% return results
res = struct('Etot',Ci(9), 'Edata',Ci(10), 'E4',Ci(11), 'E8',Ci(12),...
             'Ctot',Ci(1), 'Cdata',Ci(2), 'C4',Ci(3), 'C8',Ci(4),...
             'lambda',Ci(13), 'bF4',Ci(5), 'bF8',Ci(6), ...
             'n4',Ci(7), 'n8',Ci(8), 'dR',Ci(14), 'er',er, ...
             'minsumRi',Ci(15), 'maxsumRi',Ci(16), 'text',ti, ...
             'bF8percent',bF8percent, 'M',M, 'I',I, 'J',J, 'useMex',useMex );

if (exist('scaleRi','var')); res.scaleRi = scaleRi; end;         
if (exist('cc','var') && isstruct(cc))
    if isfield(cc, 'NumObjects'); res.ccNumObjects = cc.NumObjects; end;
    if isfield(cc, 'NumPixels'); res.ccNumPixels = cc.NumPixels; end;
    if isfield(cc, 'LenBorder'); res.ccLenBorder = cc.LenBorder; end;
    if isfield(cc, 'OutBorder'); res.ccOutBorder = cc.OutBorder; end;
end
if returnRi;  res.Ri = Ri; end;
if returnLi;  res.Li = Li; end;
if (exist('Ltrue','var') && returnLtrue); res.Ltrue = Ltrue; end;

if nargout >= 1
    varargout{1} = res;
end
if nargout >= 2
    for i=2:nargout
        varargout{i} = [];
    end
end

return

function C = LIS_cost_local(Ri, Li, p)
M=p(1); I=p(2); J=p(3); bF4=p(5); bF8=p(6); lambda=p(9);
dR=p(10); minsumRi=p(11); maxsumRi=p(12);
idxL = find(Li);   % the nonzeros in L
if (bF4 < 0)
    [temp,L] = min(reshape(Ri,M,I*J));
    minsumRi = sum(Ri( M*(idxL-1) + reshape(L(idxL),numel(idxL),1) ));
    [temp,L] = max(reshape(Ri,M,I*J));
    maxsumRi = sum(Ri( M*(idxL-1) + reshape(L(idxL),numel(idxL),1) ));
    dR = maxsumRi - minsumRi;
    bF4 = floor(0.5+lambda*dR/(I*(J-1)+(I-1)*J));
    if (bF8 > 0)
        bF8 = floor(0.5+(bF8/100)*lambda*dR/(I*(J-1)+(I-1)*J));
    else
        bF8 = 0;
    end
end
Cdata = sum(Ri( M*(idxL-1) + double(Li(idxL)) ));
L = reshape(Li, I, J);
n4 = nnz((L(:,1:(J-1))~=L(:,2:J)) & (L(:,1:(J-1))>0) & (L(:,2:J)>0)) + ...
     nnz((L(1:(I-1),:)~=L(2:I,:)) & (L(1:(I-1),:)>0) & (L(2:I,:)>0));
n8 = nnz((L(1:(I-1),1:(J-1))~=L(2:I,2:J)) & (L(1:(I-1),1:(J-1))>0) & (L(2:I,2:J)>0)) + ...
     nnz((L(2:I,1:(J-1))~=L(1:(I-1),2:J)) & (L(2:I,1:(J-1))>0) & (L(1:(I-1),2:J)>0));
C4 = n4*bF4;
C8 = n8*bF8;
Ctot = Cdata + C4 + C8;
if (dR < 1); dR=1; end;
Edata = (Cdata-minsumRi)/dR;
E4 = n4/(I*(J-1)+(I-1)*J);
if bF4 > 0
    E8 = (bF8/bF4)*n8/(I*(J-1)+(I-1)*J);  
else
    E8 = 0;
end
Etot = Edata + lambda*(E4 + E8);
C = [Ctot, Cdata, C4, C8, bF4, bF8, n4, n8, ...
     Etot, Edata, E4, E8, lambda, dR, minsumRi, maxsumRi];
disp([' bF4=',int2str(bF4),' bF8=',int2str(bF8), ...
      ' Edata=',int2str(Edata),' E4=',int2str(E4), ' E8=',int2str(E8)]);
return;

