function varargout = LIS_swap(Ri, Li, varargin)
% LIS_swap         alpha-beta-swap on label image using max-flow-algorithm 
%
%  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:
%   ----------------
%   Ls = LIS_swap(R, L);       % all options are default
%   [Ls, res] = LIS_swap(R, L, 'option',value, ...)
%   [Ls, res] = LIS_swap(R, L, Ltrue, 'option',value, ...); % including true labeling
%
%     Ls  is resulting labeling after alpha-beta-swap
%     res  is a struct with more results
%     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)
%     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)
% 
%       Options used when calculating the cost for a labeling 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
% 
%       Options controlling how many times and how alpha-beta-swap is done:
%   'maxLoop' maximum number of loops to do, each loop may do several 
%             alpha-beta-swaps. If no pixel is swapped in a loop, no more
%             loops are done.
%   'alpha'   the alpha value(s) to use in each swap. 
%   'beta'    the beta value(s) to use in each swap. 
%             Both alpha and beta may be a matrices of maxLoop columns. 
%             The non-zeros elements in each column are combined by each other, 
%             and the pairs (where beta > alpha) are processed in each loop. 
%             -1 may be given to select all combinations in random order
%   'onlyNeighbors' true (default) does only process pair where alpha and 
%             beta are neighbors, i.e. have common border.
%   'bw' or 'borderWidth' restrict the swap to be done within this
%             distance from border between alpha and beta. Default 0 means that
%             'bw' restriction is not used, all alpha and beta labels are considered.
%             It may be given as a row vector, one value for each loop.
%
%        Other options
%   'scaleRi' may be given, if not a suitable value is found if needed
%   'FigNo'   number for figure window to display results in, see LIS_plotcol.m
%             default: [0,0,0] (no plot for input or output or swap region)
%   'v' or 'verbose' may be 0/false (default), 1/true or 2 (very verbose).
% 
%   Examples:
%   ---------
%   clear all; tim = 3;  Ltrue = teximage(tim+200);
%   R = LIS_getR(tim);  [I,J,M] = size(R);  
%   sigma = 5; Lg = LIS_smooth(R, sigma);
%   sPar = struct('maxLoop',10, 'bw',[1,2,3,5,7,12,19], 'alpha',-1, 'beta',-1);
%   tic; [Ls, res] = LIS_swap(R, Lg, Ltrue, sPar, 'v',1, 'FigNo',[1,2,0] ); toc;


%----------------------------------------------------------------------
% Copyright (c) 2012.  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  31.10.2012  A version that seems to work as intended is ready
% ver. 1.1  01.11.2012  Use all pairs of labels
% Ver. 1.2  16.11.2012  changed parameters
%----------------------------------------------------------------------

% some identities that may be helpful
% i = rem(pix-1,I)+1;  j = floor((pix-1)/I)+1;  pix = (j-1)*I+i;

mfile = 'LIS_swap';

%% 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} = [];
    end
    if nargout >= 2
        varargout{2} = struct('Error',t);
    end
    return
end

%% defaults, initial values
lambda = 2;
bF4 = -1;
bF8 = -1;
bF8percent = 70;
dR = -1;
minsumRi = -1;
maxsumRi = -1;
figno = [0, 0, 0];
maxLoop = 6;
bw = [3,5,8,13,21,34];
alpha_in = -1;
beta_in = -1;
onlyNeighbors = true;
verbose = 0;
    
%% 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,'alpha')
        if isnumeric(opVal); alpha_in = opVal(1); end;
    end
    if strcmpi(opName,'beta')
        if isnumeric(opVal); beta_in = opVal(1); end;
    end
    if strcmpi(opName,'maxLoop')
        if isnumeric(opVal); maxLoop = opVal(1); end;
    end
    if strcmpi(opName,'bw') || strcmpi(opName,'borderWidth')
        if isnumeric(opVal); bw = reshape(opVal(:),1,numel(opVal)); end;
    end
    if strcmpi(opName,'onlyNeighbors') 
        if islogical(opVal); onlyNeighbors = opVal(1); end;
        if isnumeric(opVal); onlyNeighbors = (opVal(1) > 0); end;
    end
    if strcmpi(opName,'FigNo')
        if (isnumeric(opVal) && (numel(opVal) >= 3)); figno = opVal(1:3); end;
        if (isnumeric(opVal) && (numel(opVal) == 2)); figno = [opVal(1),opVal(2),0]; end;
        if (isnumeric(opVal) && (numel(opVal) == 1)); figno = [0,opVal(1),0]; end;
    end
    if strcmpi(opName,'verbose') || strcmpi(opName,'v')
        if (islogical(opVal) && opVal); verbose = 1; end;
        if isnumeric(opVal); verbose = opVal(1); 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)
        [~,Li] = min(reshape(Ri, M, I*J));   % integers in Ri may give a little
        % 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 some important parameters
if (maxLoop < 1); maxLoop = 1; end;
n = min(maxLoop,numel(bw));
bw = [reshape(bw(1:n),1,n), zeros(1,maxLoop-n)];
if (alpha_in(1) < 0)
    alpha_in = zeros(M,maxLoop);
    for i=1:maxLoop
        alpha_in(:,i) = reshape(randperm(M),M,1);
    end
end
if (beta_in(1) < 0)
    beta_in = zeros(M,maxLoop);
    for i=1:maxLoop
        beta_in(:,i) = reshape(randperm(M),M,1);
    end
end

%% Display info for input
if (verbose > 1)  % very verbose
    disp(' ');
    disp([mfile,' started ',datestr(now),'  (lambda=',num2str(lambda),')']);
    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)),')']);
    end
    if (maxLoop > 1)
        disp([' Do up to ',int2str(maxLoop),' loops through all labels.']);
        if nnz(bw)
            disp([' (bw: ',nnz(bw),' non-zeros from ',num2str(min(nonzeros(bw))), ...
                ' to ',num2str(max(bw)),')']);
        end
    elseif (numel(alpha_in) == 1) && (numel(beta_in) == 1) 
        disp([' Do one alpha-beta-swap (alpha=',int2str(alpha_in),...
            ', beta=',int2str(beta_in),', bw=',num2str(bw(1)),').']);
    end
end
% disp([' ------ verbose = ',int2str(verbose)]);

%% do the calulations
paramCost = struct('M',M, 'I',I, 'J',J, 'bF4',bF4, 'bF8',bF8,...
    'dR',dR, 'minsumRi',minsumRi, 'maxsumRi',maxsumRi, 'verbose',0, 'FigNo',0, ...
    'bF8percent',bF8percent, 'lambda',lambda, 'useMex',true);
if exist('Ltrue','var') 
    costLtrue = LIS_cost(Ri, Ltrue, Ltrue, paramCost);      % cost of input labeling
    costLtrue.text(10:17) = '   Ltrue';
    paramCost.bF4 = costLtrue.bF4; 
    paramCost.bF8 = costLtrue.bF8; 
    paramCost.lambda = costLtrue.lambda; 
    paramCost.dR = costLtrue.dR; 
    paramCost.minsumRi = costLtrue.minsumRi; 
    paramCost.maxsumRi = costLtrue.maxsumRi;
    paramCost.FigNo = figno(1);
    costLi = LIS_cost(Ri, Li, Ltrue, paramCost);      % cost of input labeling
else
    paramCost.FigNo = figno(1);
    costLi = LIS_cost(Ri, Li, paramCost);      % cost of input labeling
    paramCost.bF4 = costLi.bF4; 
    paramCost.bF8 = costLi.bF8; 
    paramCost.lambda = costLi.lambda;
    paramCost.dR = costLi.dR; 
    paramCost.minsumRi = costLi.minsumRi; 
    paramCost.maxsumRi = costLi.maxsumRi;
end
bF4 = costLi.bF4;
bF8 = costLi.bF8;
lambda = costLi.lambda;
if (verbose >= 1)
    disp([' The used values for bF4 and bF8 and lambda are:  bF4=',num2str(bF4),...
            ', bF8=',num2str(bF8),' and lambda=',num2str(lambda),'.']);
    if exist('Ltrue','var') 
        disp([' ',costLtrue.text]);    
    end
    disp([' ',costLi.text]);    
end

L = reshape(Li,I,J);   % IxJ matrix of uint8
n4hood = strel([0,1,0;1,1,1;0,1,0]);  


% alpha-beta-swap for given pairs of labels 
for loop = 1:maxLoop
    totalSwapped = 0;
    if (bw(loop) > 0)       
        nh1 =  strel('disk',floor(bw(loop)),0);
    end
    for alpha = reshape(alpha_in(:,loop),1,size(alpha_in,1));
        if ((alpha < 1) || (alpha > M))
            continue; 
        end
        for beta = reshape(beta_in(:,loop),1,size(beta_in,1));
            if ((beta <= alpha) || (beta > M))
                continue; 
            end
            L1 = (L == alpha);
            L2 = (L == beta);
            if onlyNeighbors  
                commonEdges = nnz(imdilate(L1,n4hood) & L2);
                if (commonEdges == 0)
                    % disp('no common edges');
                    continue;
                end
            end
            if (bw(loop) <= 0) || (bw(loop) > max(I,J))     
                mask = (L1 | L2);
            else
                % do only consider pixels near the edge between these two labels
                mask = ((imdilate(L1,n4hood) & L2) | (imdilate(L2,n4hood) & L1)); % border pixels
                mask = imdilate(mask, nh1) & (L1 | L2);   % dilated  
            end
            
            if (verbose >= 2)
                disp([mfile,': Process ',int2str(nnz(mask)),' pixels of label ',...
                    int2str(alpha),' and ',int2str(beta)]);
            end
            
            [A, T] = abSwap_makeGraph(Ri, L, mask, bF4, bF8, alpha, beta);

            [~, labels] = LIS_maxflow(A, T);
            
            pix = find(mask);
            nofSwapped = nnz( (L(pix)==beta) ~= labels );
            totalSwapped = totalSwapped + nofSwapped;

            L(pix(labels == 0)) = alpha;
            L(pix(labels == 1)) = beta;
            
            if (verbose >= 2)
                disp(['  number of pixels swapped is ',int2str(nofSwapped)]);
            end
        end
    end
    if (verbose >= 1)
        disp(['  Total number of pixels swapped in loop ',int2str(loop),...
            '    is    ',int2str(totalSwapped)]);
    end
    if (totalSwapped == 0)
        break;
    end
end

if ((figno(3) > 0) && exist('mask', 'var'))
    figure(figno(3)); spy(mask);
    title([strrep(mfile,'_','-'),': Pixels processed in loop ',int2str(loop),...
        ', \alpha=',int2str(alpha),', \beta=',int2str(beta)]);
end

paramCost.FigNo = figno(2);
if exist('Ltrue','var') 
    costLs = LIS_cost(Ri, L(:), Ltrue, paramCost);
else
    costLs = LIS_cost(Ri, L(:), paramCost);  
end
costLs.text(10:17) = '   Lswap';
if (verbose >= 1)
    disp([' ',costLs.text]);    
end


%% return results
res = struct('costLi',costLi, 'costLs',costLs, 'paramCost',paramCost);

if (exist('costTrue','var')); res.costTrue = costTrue; end;         

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

return


%%  make the graph
function [A, T] = abSwap_makeGraph(Ri, L, mask, bF4, bF8, alpha, beta)
% build the graph, edges between nodes in A, and in T from s and to t
[I,J] = size(L);
M = numel(Ri)/(I*J);
nnzPix = nnz(mask);
pix = find(mask);  % note pix here is not as in alphaExpand_makeGraph

pixL = pix-I;  pixL(pix <= I) = 0;  % pixL = neighbor_left(pix);
pixA = pix-1;  pixA(rem(pix,I) == 1) = 0;  % pixA = neighbor_above(pix);
if (bF8 > 0)
    pixLA = pix-I-1;
    pixLA((pix <= I)|(rem(pix,I) == 1)) = 0;
    pixLB = pix-I+1;
    pixLB((pix <= I)|(rem(pix,I) == 0)) = 0;
end
% pixNode = reshape(1:nnzPix,nnzPix,1);
idx = zeros(I,J);   % idx is like pix in alphaExpand_makeGraph
idx(pix) = reshape(1:nnzPix,nnzPix,1);
pixNodeL = zeros(nnzPix,1);
pixNodeL(pixL>0) = idx(nonzeros(pixL));
pixNodeA = zeros(nnzPix,1); pixNodeA(pixA>0) = idx(nonzeros(pixA));
nz4 = nnz(pixNodeL)+nnz(pixNodeA);
if (bF8 > 0)
    pixNodeLA = zeros(nnzPix,1);
    pixNodeLA(pixLA>0) = idx(nonzeros(pixLA));
    pixNodeLB = zeros(nnzPix,1); pixNodeLB(pixLB>0) = idx(nonzeros(pixLB));
    nz8 = nnz(pixNodeLA)+nnz(pixNodeLB);
    i = [nonzeros(pixNodeL);nonzeros(pixNodeA);nonzeros(pixNodeLA);nonzeros(pixNodeLB)];
    j = [find(pixNodeL);find(pixNodeA);find(pixNodeLA);find(pixNodeLB)];
    A = sparse([i;j],[j;i],...
        [bF4*ones(nz4,1);bF8*ones(nz8,1);bF4*ones(nz4,1);bF8*ones(nz8,1)],...
        nnzPix,nnzPix,2*(nz4+nz8));
else
    i = [nonzeros(pixNodeL);nonzeros(pixNodeA)];
    j = [find(pixNodeL);find(pixNodeA)];
    A = sparse([i;j], [j;i], bF4*ones(2*nz4,1), nnzPix, nnzPix, 2*nz4);
end

i = M*(pix-1)+double(alpha);
j = M*(pix-1)+double(beta);
T2 = double(Ri(i));
T1 = double(Ri(j));
% should adjust T for pixels/nodes near the edges where the outside is alpha or beta
outside = find(xor(imdilate(mask, strel([1,1,1;1,1,1;1,1,1])),mask));
p = outside((L(outside) == alpha)); 
% if an outside neighbor pixel is alpha, we should be less willing to label the
% border pixel as betaa and thus increase cost for selecting beta, Ri(j), which
% here is in T2.   ????????????????
if numel(p) > 0
    pn = [nonzeros(max(p-I,0)); ...                   % left
          nonzeros((p-1).*double(rem(p,I)~=1)); ...   % above
          nonzeros((p+1).*double(rem(p,I)~=0)); ...   % below
          nonzeros((p+I).*double(p<=(I*J-I))) ];      % right
    nodes = idx(pn(mask(pn))); % the nodes that have outside neighbors equal alpha
    for i=1:numel(nodes); T1(nodes(i))=T1(nodes(i))+bF4; end;
    if (bF8 > 0)
        pn = [nonzeros(max(p-I-1,0).*double(rem(p,I)~=1)); ...   % left above
              nonzeros(max(p-I+1,0).*double(rem(p,I)~=0)); ...   % left below
              nonzeros((p+I-1).*double((rem(p,I)~=1)&(p<=(I*J-I)))); ...   % right above
              nonzeros((p+I+1).*double((rem(p,I)~=0)&(p<=(I*J-I))))];   % right below
        nodes = idx(pn(mask(pn))); % the nodes that have outside neighbors equal alpha
        for i=1:numel(nodes); T1(nodes(i))=T1(nodes(i))+bF8; end;
    end
end
p = outside((L(outside) == beta)); 
if numel(p) > 0
    pn = [nonzeros(max(p-I,0)); ...                   % left
          nonzeros((p-1).*double(rem(p,I)~=1)); ...   % above
          nonzeros((p+1).*double(rem(p,I)~=0)); ...   % below
          nonzeros((p+I).*double(p<=(I*J-I))) ];      % right
    nodes = idx(pn(mask(pn))); % the nodes that have outside neighbors equal beta
    for i=1:numel(nodes); T2(nodes(i))=T2(nodes(i))+bF4; end;
    if (bF8 > 0)
        pn = [nonzeros(max(p-I-1,0).*double(rem(p,I)~=1)); ...   % left above
              nonzeros(max(p-I+1,0).*double(rem(p,I)~=0)); ...   % left below
              nonzeros((p+I-1).*double((rem(p,I)~=1)&(p<=(I*J-I)))); ...   % right above
              nonzeros((p+I+1).*double((rem(p,I)~=0)&(p<=(I*J-I))))];   % right below
        nodes = idx(pn(mask(pn))); % the nodes that have outside neighbors equal beta
        for i=1:numel(nodes); T2(nodes(i))=T2(nodes(i))+bF8; end;
    end
end

% cancel direct flow from S to pixel/node to T
T = sparse([max(T1-T2,0),max(T2-T1,0)]);

return
