function varargout = EMP8_swap(s)
% EMP8_swap        alpha-beta-swap on label image using max-flow-algorithm 
%                  EMP8: Energy Minimization using Potts model and 8 neighbors
%
% This function may use:
%   LIS_maxflow(..)       call maxflowmex (a maxflow implementation of 
%                         Boykov-Kolmogorov's algorithm by Michael Rubinstein)
%   EMP8_cost_mex(..)     to find cost for a labeling
%
% Use of function:
% ----------------
%   L = EMP8_swap(s)
%   [L,c] = EMP8_swap(s)
% 
%     L   is the resulting new labeling after swap moves
%     c   is a struct with cost results in each field
%     s   is a struct as returned by EMP8_argin, see help for that file
%    Special options, i.e. fields in s 
%     s.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.
%     s.alpha  the alpha value(s) to use in each swap. 
%     s.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 if 
%         s.useRand is true, else just in 'sorted' order
%     s.onlyNeighbors   true (default) does only process pairs where alpha and 
%         beta are neighbors, i.e. have common border.
%     s.borderWidth   restrict the swap to be done within this
%         distance from border between alpha and beta. Default 0 means that
%         this restriction is not used, all alpha and beta labels are considered.
%         It may be given as a row vector, one value for each loop.
% 
% Example:
% -------------------------------------------------------
%  clear all; tim = 4;  sigma = 5;  
%  Uim = teximage(100+tim,[],0,4);  % read test image, and display it
% 
%  s = EMP8_argin(tim, sigma, tim, 'v',2, 'bw',-1, 'useRand',false); 
%  [L,c] = EMP8_swap(s);
%  figure(1); clf;  % and plot results in figure 1
%  figTitle = ['Labeling of test image ',int2str(tim),' after using EMP8-swap, swUpq=1.'];
%  LIS_plotcol(reshape(L,s.I,s.J), figTitle, c.errorRate, -1, c.Etot, c.Edata, c.E4+c.E8);
%
%  [I,J] = size(Uim); Upq = zeros(I,J,4);
%  Upq(:,:,1) = Uim - [Uim(1,:); Uim(2:I,1),Uim(1:(I-1),1:(J-1))];
%  Upq(:,:,2) = Uim - [Uim(:,1),Uim(:,1:(J-1))];
%  Upq(:,:,3) = Uim - [Uim(:,1),[Uim(2:I,1:(J-1));Uim(I,2:J)]];
%  Upq(:,:,4) = Uim - [Uim(1,:);Uim(1:(I-1),:)];
%  Upq = abs(Upq); Upq = 1 - Upq/max(Upq(:));
%  s = EMP8_argin(tim, sigma, tim, 'swUpq',2, 'Upq',Upq, 'v',2, 'useRand',false);
%  [L,c] = EMP8_swap(s);
%  figure(2); clf;  % and plot results in figure 2
%  figTitle = ['Labeling of test image ',int2str(tim),' after using EMP8-swap, swUpq=2.'];
%  LIS_plotcol(reshape(L,s.I,s.J), figTitle, c.errorRate, -1, c.Etot, c.Edata, c.E4+c.E8);
%
%  s = EMP8_argin(tim, sigma, tim, 'swUpq',3, 'Uim',tim, 'Utab',[2,120,0.2], 'v',2, 'bw',-1, 'useRand',false);
%  [L,c] = EMP8_swap(s);
%  figure(3); clf; 
%  figTitle = ['Labeling of test image ',int2str(tim),' after using EMP8-swap, swUpq=3.'];
%  LIS_plotcol(reshape(L,s.I,s.J), figTitle, c.errorRate, -1, c.Etot, c.Edata, c.E4+c.E8);


%----------------------------------------------------------------------
% 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
% Ver. 1.3  12.12.2013  EMP8_swap made based on LIS_swap
% Ver. 1.4  14.10.2014  some minor modifications
%----------------------------------------------------------------------

if ( ~exist('EMP8_cost_mex','file') )
    c = [mfilename,': 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
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.bF4 < 0) || (s.bF8 < 0))
    % this should normally not happen
    disp([mfilename,': 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([mfilename,': 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([mfilename,': some values are bF4 = ',int2str(s.bF4),', bF8 = ',int2str(s.bF8),...
          ', lambda = ',num2str(s.lambda),', bF8percent = ',int2str(s.bF8percent), ...
          ', swUpq = ',int2str(s.swUpq)]);
end

%% check some important parameters
I = s.I; 
J = s.J;
M = s.M;
if (s.borderWidth(1) == -1); 
    bw = [3,5,8,13,21,34]; 
else
    bw = s.borderWidth;
end
maxLoop = s.maxLoop;
if (maxLoop < 1); maxLoop = 1; end;
n = min(maxLoop,numel(bw));
bw = [reshape(bw(1:n),1,n), zeros(1,maxLoop-n)];
if (s.alpha(1) < 0)
    if s.useRand
        alpha_in = zeros(M,maxLoop);
        for i=1:maxLoop
            alpha_in(:,i) = reshape(randperm(M),M,1);
        end
    else
        alpha_in = repmat(reshape(1:M,M,1),1,maxLoop);
    end
else
    alpha_in = s.alpha;
end
if (s.beta(1) < 0)
    if s.useRand
        beta_in = zeros(M,maxLoop);
        for i=1:maxLoop
            beta_in(:,i) = reshape(randperm(M),M,1);
        end
    else
        beta_in = repmat(reshape(1:M,M,1),1,maxLoop);
    end
else
    beta_in = s.beta;
end

%% the actual work her is to prepare for maxflow algorithm, i.e. build the graph
tic;
L = reshape(s.Li,I,J);
n4hood = strel([0,1,0;1,1,1;0,1,0]);  

% alpha-beta-swap for given pairs of labels 
for loop = 1:maxLoop
    if (s.verbose == 2)
       disp([mfilename,': loop ',int2str(loop),', bw(',int2str(loop),') = ',num2str(bw(loop)),'.']);
    end
    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 s.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 (s.verbose >= 3)
                disp([mfilename,': Process ',int2str(nnz(mask)),' pixels of label ',...
                    int2str(alpha),' and ',int2str(beta)]);
            end
            
            [A, T] = abSwap_makeGraph(s, L, mask, 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 (s.verbose == 2)
                disp([' processed ',sprintf('%7i',nnz(mask)),' pixels of label ',...
                    sprintf('%2i',alpha),' and ',sprintf('%2i',beta),' and swapped ',...
                    sprintf('%7i',nofSwapped),' pixels.']);
            end
            if (s.verbose >= 3)
                disp(['  number of pixels swapped is ',int2str(nofSwapped)]);
            end
        end
    end
    if (s.verbose >= 1)
        disp([' Total number of pixels swapped in loop ',sprintf('%2i',loop),...
            ' is  ',sprintf('%7i',totalSwapped)]);
    end
    if (totalSwapped == 0)
        break;
    end
end
L = L(:);
swapTime = toc;


%% return results
if (nargout >= 1)
    varargout{1} = L;  
end
if (nargout >= 2)
    updated = 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
    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
    swapCost = sprintf(s.formatCost, 'Lswap', c([10,13,11,12,9]), 100*er);
    if (er < 0) % it is nicer to print nothing rather the negative error rate
        swapCost = [swapCost(1:72),'.'];
    end
    if (s.verbose >= 1)
        disp([mfilename,': done. Swapped total ',int2str(updated),', pixels using ',...
              num2str(swapTime),' seconds.']);
        disp(swapCost);
    end
    % den normale situasjonen
    c(17:20) = [updated; swapTime; -1; er];
    c = costTilStruct(c);
    c.swapCost = swapCost; 

    varargout{2} = c;
end

return

%% funksjoner

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),'unused',c(19),'errorRate',c(20));
return;

function [A, T] = abSwap_makeGraph(s, L, mask, alpha, beta)
% build the graph, edges between nodes in A, and in T from s and to t
I = s.I; 
J = s.J;
M = s.M;
% [I,J] = size(L);    % mask is also size IxJ (logical)
% M = numel(Ri)/(I*J);
nnzPix = nnz(mask);
pix = find(mask);     % a column vector

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 (s.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(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 (s.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);
    if (s.swUpq == 1)
        i = [nonzeros(pixNodeL);nonzeros(pixNodeA);nonzeros(pixNodeLA);nonzeros(pixNodeLB)];
        j = [find(pixNodeL);find(pixNodeA);find(pixNodeLA);find(pixNodeLB)];
        Aval = [s.bF4*ones(nz4,1);s.bF8*ones(nz8,1)];
    elseif (s.swUpq == 2)
        i = [nonzeros(pixNodeL);nonzeros(pixNodeA);nonzeros(pixNodeLA);nonzeros(pixNodeLB)];
        jL = find(pixNodeL);
        jA = find(pixNodeA);
        jLA = find(pixNodeLA);
        jLB = find(pixNodeLB);
        j = [jL; jA; jLA; jLB];
        Aval = [(s.bF4/s.PU1) * double(s.Upq( 4*pix(jL)-2 )); ...
                (s.bF4/s.PU1) * double(s.Upq( 4*pix(jA) )); ...
                (s.bF8/s.PU1) * double(s.Upq( 4*pix(jLA)-3 )); ...
                (s.bF8/s.PU1) * double(s.Upq( 4*pix(jLB)-1 ))   ]; 
    else % (s.swUpq == 3)
        i4 = [nonzeros(pixNodeL);nonzeros(pixNodeA)];
        j4 = [find(pixNodeL);find(pixNodeA)];
        i8 = [nonzeros(pixNodeLA);nonzeros(pixNodeLB)];
        j8 = [find(pixNodeLA);find(pixNodeLB)];
        i = [i4; i8];
        j = [j4; j8];
        Aval = [(s.bF4/s.PU1) * double(s.Utab( abs(double(s.Uim(pix(i4)))-double(s.Uim(pix(j4))))+1 )); ...
                (s.bF8/s.PU1) * double(s.Utab( abs(double(s.Uim(pix(i8)))-double(s.Uim(pix(j8))))+1 )) ]; 
    end
    A = sparse([i;j], [j;i], [Aval;Aval], nnzPix, nnzPix, 2*(nz4+nz8));
else
    i = [nonzeros(pixNodeL);nonzeros(pixNodeA)];
    if (s.swUpq == 1)
        j = [find(pixNodeL);find(pixNodeA)];
        Aval = s.bF4*ones(nz4,1);
    elseif (s.swUpq == 2)
        jL = find(pixNodeL);
        jA = find(pixNodeA);
        j = [jL; jA];
        Aval = [(s.bF4/s.PU1) * double(s.Upq( 4*pix(jL)-2 )); ...
                (s.bF4/s.PU1) * double(s.Upq( 4*pix(jA) )) ]; 
    else % (s.swUpq == 3)
        j = [find(pixNodeL);find(pixNodeA)];
        Aval = (s.bF4/s.PU1) * double(s.Utab( abs(double(s.Uim(pix(i)))-double(s.Uim(pix(j))))+1 ));
    end
    A = sparse([i;j], [j;i], [Aval;Aval], nnzPix, nnzPix, 2*nz4);
end

i = M*(pix-1)+double(alpha);
j = M*(pix-1)+double(beta);
T2 = double(s.Ri(i));
T1 = double(s.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));

for ab = [alpha, beta]
    % disp(['ab = ',int2str(ab),' (',class(ab),')']); 
    p = outside((L(outside) == ab));
    % if an outside neighbor pixel is alpha, we should be less willing to label the
    % border pixel as beta and thus increase cost for selecting beta, Ri(j),
    % this cost is now in T2.
    if numel(p) > 0
        if (s.swUpq == 1)
            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
            if (numel(nodes) > 0)
                if (ab == alpha)
                    T1(nodes) = T1(nodes) + s.bF4;
                else % (ab == beta)
                    T2(nodes) = T2(nodes) + s.bF4;
                end
            end
            if (s.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
                if (numel(nodes) > 0)
                    if (ab == alpha)
                        T1(nodes) = T1(nodes) + s.bF8;
                    else  % (ab == beta)
                        T2(nodes) = T2(nodes) + s.bF8;
                    end
                end
            end
        elseif (s.swUpq == 2)  % use s.Upq
            pn = max(p-I,0);             % left of p (p in (outside == ab))
            p2 = p(pn > 0);              % p where pn exist
            pn = nonzeros(pn);           % and these pn values
            % mask(pn) indicate which values of pn and p2 (or p) that are relevant
            nodes = idx(pn(mask(pn)));   % nodes (in T1/T2) for relevant pn
            if (numel(nodes) > 0)
                p2 = p2(mask(pn));       % and index(es) for relevant p (p2)
                if (ab == alpha)
                    T1(nodes) = T1(nodes) + (s.bF4/s.PU1) * double(s.Upq( 4*p2-2 ));
                else % (ab == beta)
                    T2(nodes) = T2(nodes) + (s.bF4/s.PU1) * double(s.Upq( 4*p2-2 ));
                end
            end
            % and above
            pn = (p-1).*double(rem(p,I)~=1);
            p2 = p(pn > 0);
            pn = nonzeros(pn);
            nodes = idx(pn(mask(pn)));
            if (numel(nodes) > 0)
                p2 = p2(mask(pn));     
                if (ab == alpha)
                    T1(nodes) = T1(nodes) + (s.bF4/s.PU1) * double(s.Upq( 4*p2 ));
                else % (ab == beta)
                    T2(nodes) = T2(nodes) + (s.bF4/s.PU1) * double(s.Upq( 4*p2 ));
                end
            end
            % and below
            pn = nonzeros((p+1).*double(rem(p,I)~=0));
            p2 = pn(mask(pn));    % here p2 is not from p but from pn!
            if (numel(p2) > 0)
                nodes = idx(p2); 
                if (ab == alpha)
                    T1(nodes) = T1(nodes) + (s.bF4/s.PU1) * double(s.Upq( 4*p2 ));
                else % (ab == beta)
                    T2(nodes) = T2(nodes) + (s.bF4/s.PU1) * double(s.Upq( 4*p2 ));
                end                
            end
            % and right
            pn = nonzeros((p+I).*double(p<=(I*J-I)));
            p2 = pn(mask(pn));    % here p2 is not from p but from pn!
            if (numel(p2) > 0)
                nodes = idx(p2); 
                if (ab == alpha)
                    T1(nodes) = T1(nodes) + (s.bF4/s.PU1) * double(s.Upq( 4*p2-2 ));
                else % (ab == beta)
                    T2(nodes) = T2(nodes) + (s.bF4/s.PU1) * double(s.Upq( 4*p2-2 ));
                end                
            end
            if (s.bF8 > 0)
                pn = max(p-I-1,0).*double(rem(p,I)~=1);  % left above
                p2 = p(pn > 0);
                pn = nonzeros(pn);
                nodes = idx(pn(mask(pn)));
                if (numel(nodes) > 0)
                    p2 = p2(mask(pn));     
                    if (ab == alpha)
                        T1(nodes) = T1(nodes) + (s.bF8/s.PU1) * double(s.Upq( 4*p2-3 ));
                    else % (ab == beta)
                        T2(nodes) = T2(nodes) + (s.bF8/s.PU1) * double(s.Upq( 4*p2-3 ));
                    end
                end
                % and left below
                pn = max(p-I+1,0).*double(rem(p,I)~=0);
                p2 = p(pn > 0);
                pn = nonzeros(pn);
                nodes = idx(pn(mask(pn)));
                if (numel(nodes) > 0)
                    p2 = p2(mask(pn));     
                    if (ab == alpha)
                        T1(nodes) = T1(nodes) + (s.bF8/s.PU1) * double(s.Upq( 4*p2-1 ));
                    else % (ab == beta)
                        T2(nodes) = T2(nodes) + (s.bF8/s.PU1) * double(s.Upq( 4*p2-1 ));
                    end
                end
                % and right above
                pn = nonzeros((p+I-1).*double((rem(p,I)~=1)&(p<=(I*J-I))));
                p2 = pn(mask(pn));    % here p2 is not from p but from pn!
                if (numel(p2) > 0)
                    nodes = idx(p2); 
                    if (ab == alpha)
                        T1(nodes) = T1(nodes) + (s.bF8/s.PU1) * double(s.Upq( 4*p2-1 ));
                    else % (ab == beta)
                        T2(nodes) = T2(nodes) + (s.bF8/s.PU1) * double(s.Upq( 4*p2-1 ));
                    end                
                end
                % and right below
                pn = nonzeros((p+I+1).*double((rem(p,I)~=0)&(p<=(I*J-I))));
                p2 = pn(mask(pn));    % here p2 is not from p but from pn!
                if (numel(p2) > 0)
                    nodes = idx(p2); 
                    if (ab == alpha)
                        T1(nodes) = T1(nodes) + (s.bF8/s.PU1) * double(s.Upq( 4*p2-3 ));
                    else % (ab == beta)
                        T2(nodes) = T2(nodes) + (s.bF8/s.PU1) * double(s.Upq( 4*p2-3 ));
                    end                
                end
            end
        else  % (s.swUpq == 3)  use s.Uim and s.Utab
            pn = [max(p-I,0); ...                   % left
                  (p-1).*double(rem(p,I)~=1); ...   % above
                  (p+1).*double(rem(p,I)~=0); ...   % below
                  (p+I).*double(p<=(I*J-I)) ];      % right
            p2 = repmat(p,4,1);  
            p2 = p2(pn > 0);
            pn = nonzeros(pn);
            nodes = idx(pn(mask(pn))); % the nodes that have outside neighbors equal alpha
            p2 = p2(mask(pn));
            if (numel(nodes) > 0)
                % dispVar('T1   ',T1); dispVar('pix  ',pix);    
                % dispVar('nodes',nodes); dispVar('p2   ',p2);
                if (ab == alpha)
                    T1(nodes) = T1(nodes) + (s.bF4/s.PU1) * ...
                        double(s.Utab( abs(double(s.Uim(pix(nodes)))-double(s.Uim(p2)))+1 ));
                else % (ab == beta)
                    T2(nodes) = T2(nodes) + (s.bF4/s.PU1) * ...
                        double(s.Utab( abs(double(s.Uim(pix(nodes)))-double(s.Uim(p2)))+1 ));
                end
            end
            if (s.bF8 > 0)
                pn = [max(p-I-1,0).*double(rem(p,I)~=1); ...   % left above
                      max(p-I+1,0).*double(rem(p,I)~=0); ...   % left below
                      (p+I-1).*double((rem(p,I)~=1)&(p<=(I*J-I))); ...   % right above
                      (p+I+1).*double((rem(p,I)~=0)&(p<=(I*J-I)))];   % right below
                p2 = repmat(p,4,1);  
                p2 = p2(pn > 0);
                pn = nonzeros(pn);
                nodes = idx(pn(mask(pn))); % the nodes that have outside neighbors equal alpha
                p2 = p2(mask(pn));
                if (numel(nodes) > 0)
                    if (ab == alpha)
                        T1(nodes) = T1(nodes) + (s.bF8/s.PU1) * ...
                            double(s.Utab( abs(double(s.Uim(pix(nodes)))-double(s.Uim(p2)))+1 ));
                    else  % (ab == beta)
                        T2(nodes) = T2(nodes) + (s.bF8/s.PU1) * ...
                            double(s.Utab( abs(double(s.Uim(pix(nodes)))-double(s.Uim(p2)))+1 ));
                    end
                end
            end
        end
    end
end

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

return

% function dispVar(name,x)
% disp([name,' [numel, min, max] = ',int2str([numel(x), min(x), max(x)])]);
% return
    
