% Implementation and Complexity of the
% Watershed-from-Markers Algorithm Computed as a Minimal
% Cost Forest
% Petr Felkel
%
function [c] = ift(varargin)

u = varargin{1};
minima = varargin{2};
% GRADTYPE = 1: Normal gradient
% GRADTYPE = 2: Lower slope derivations as in 
%
% The Watershed Transform: Definitions, Algorithms and
% Jos B.T.M. Roerdink and Arnold Meijster
% Parallelization Strategies
% Fundamenta Informaticae 41 (2001) 187{228 1
% IOS Press
gradtype = 1;
conn = 4;
for i = 3 : 2 : nargin
    namehere = varargin{i};
    varhere = varargin{i+1};
    switch(namehere)
        case('conn')
            conn = varhere;
        case('gradtype')
            gradtype = varhere;
            if ne(sum(gradtype == [1 2]),1);
                error('Wrong value of GRADTYPE')
            end;
        otherwise
            error('Wrong input')
    end;
end;


if isempty(find(minima))
    error('No minima given')
end;

% u = imfilter(u,fspecial('gaussian',5,3),'replicate');

% compute lower slope if needed
if gradtype == 2
    disp('Computing lower slope')
    lowsl = lowerslope(u,conn);
end;

disp('Computing topograhpical distance function');

[M N] = size(u);

% intilate all pixels in the flag array(!), not image
% flagarray = (1:M*N)';
flag = zeros(M,N);

% Initialize the cost array, 
% all non-marker nodes are infinity, the others are zero
c = Inf(M,N);
ind = minima == 1;
c(ind) = u(ind);

% put the markers in the array queue
queueind = find(minima == 1);
queueval = u(minima == 1);

while 1   
    
    % get the minimum value from the queue
    [cv,indminval] = min(queueval);
    
    %
    % Find the pixel with the minimum value of the queue and take it out
    %

    % label this pixel as DONE and remove it from the queue
    indv = queueind(indminval);
    flag(indv) = 1;
    queueind(indminval) = [];
    queueval(indminval) = [];
    
    % find the neighbours of this pixel inside flag array, therefore we go to
    % coordinates representation
    [x y] = ind2sub([M N],indv);    
    
    % 8 neighbourhood
    if conn == 4        
        % 4 neighbourhood
        xp = [x x-1 x+1 x]';
        yp = [y-1 y y y+1]';
    elseif conn == 8
        xp = [x-1 x x+1 x-1 x+1 x-1 x x+1]';
        yp = [y-1 y-1 y-1 y y y+1 y+1 y+1]';
    end;

    % remove out of range subscripts
    ind = find(xp < 1 | xp > M | yp < 1 | yp > N);
    xp(ind) = [];
    yp(ind) = [];
    indp = sub2ind([M N],xp,yp);
    
    % only the pixels that are inside the FLAG array        
    insideflag = flag(indp);
    indp = indp(insideflag == 0);
    
    % we must update the coordinates also since we may have removed some
    % points that were not in FLAGARRAY
    [xp yp] = ind2sub([M N],indp);
    
    % loop over the neighbours of the current pixel that are inside the flag array
    numindp = length(indp);

    % the Euclidean distance from v to p pixels to use in the gradient
    % computation
    dist = sqrt((x*ones(numindp,1)-xp).^2 + (y*ones(numindp,1)-yp).^2);

    for i = 1 : numindp
        % the indices of the current neighbour
        indphere = indp(i);
        
        % the distance from the center pixel to thie neighbour
        disthere = dist(i);
        
        % the cost of the current neighbour
        cphere = c(indphere);

        % get gradient
        if gradtype == 1            
            cvp = getgrad(u,indphere,indv,disthere);
        elseif gradtype == 2
            cvp = getgradls(u,indphere,indv,disthere,lowsl);
        else
            error('Wrong gradtype')
        end;

        % test if its che   aper to walk to p via v compared to the old value
        costviav = cv + cvp;
        
        % if the cost via v tp p is smaller than the already given cost of
        % p then update to the lower cost
        if costviav < cphere
            
            % update the cost function
            c(indphere) = costviav;
            
            % find if p was already in the queue
            ind = find(queueind == indphere);

            % if p is in the queue from before than remove it from it
            if ~isempty(ind)
                % remove from queue
                queueind(ind) = [];
                queueval(ind) = [];
            end;
                            
            % add p back into the queue
            queueind = [queueind ; indphere];
            queueval = [queueval ; costviav];            
        end;

%         uhere = u;
%         uhere(indv) = min(u(:));
%         uhere(indphere) = max(u(:));
%         uhere2 = u;
%         uhere2(queueind) = max(u(:));
%         show(uhere,1)
%         show(flag,2)
%         show(c,3)
%         show(u,4)
%         show(uhere2,5)
%         numindp
%         i
%         indv
%         indphere
%         ~isempty(ind)
%         queueind
%         queueval        
%         costviav
%         cphere
%         pause
%         
    end;
%     queueval
%     uhere2 = u;
%     uhere2(queueind) = max(u(:));   
%     show(c,3)
%     show(u,4)
%     show(uhere2,5)
%     pause
    
    % stop if queue is empty
    if isempty(queueind)
        break;
    end;
    
end;

% show(u,1)
% show(c,3)

%------------------------------------------------

function [cvp] = getgrad(u,indp,indv,dist)


% p : neighbour
% v : this pixel

% the gradient, slope
cvp = abs(u(indp) - u(indv))/dist;

% ----------------------------------------------

function [cvp] = getgradls(u,indp,indv,dist,lowsl)
           


% v : this pixel
uv = u(indv);

% p : neighbour
up = u(indp);

if uv > up
    % lower slope of this pixel
    cvp = lowsl(indv)*dist;
elseif uv < up
    % lower slope of neighbour pixel
    cvp = lowsl(indp)*dist;       
elseif up == uv
    % combination
    cvp = (1/2)*(lowsl(indv) + lowsl(indp))*dist;
end;        

