% CLEANEDGELIST - remove short edges from a set of edgelists
%
% Function to clean up a set of edge lists generated by EDGELINK so that
% isolated edges and spurs that are shorter that a minimum length are removed.
% This code can also be use with a set of line segments generated by LINESEG.
%
% Usage: nedgelist = cleanedgelist(edgelist, minlength)
%
% Arguments:
%                edgelist - a cell array of edge lists in row,column coords in
%                           the form
%                           { [r1 c1   [r1 c1   etc }
%                              r2 c2    ...
%                                       ...
%                              rN cN]   ....]   
%                minlength - minimum edge length of interest
%
% Returns:
%                nedgelist - the new, cleaned up set of edgelists
%
% See also:  EDGELINK, DRAWEDGELIST, LINESEG

% Copyright (c) 2006-2007 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% pk at csse uwa edu au
% http://www.csse.uwa.edu.au/
% 
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
% 
% The above copyright notice and this permission notice shall be included in 
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.
%
% December 2006  Original version
% February 2007  A major rework to fix several problems, hope they really are
%                fixed!


function nedgelist = cleanedgelist(edgelist, minlength)
    
    Nedges = length(edgelist);
    Nnodes = 2*Nedges;

    % Each edgelist has two end nodes - the starting point and the ending point.
    % We build up an adjacency/connection matrix for each node so that we can
    % determine which, if any, edgelists are connected to a node. We also
    % maintain an adjacency matrix for the edges themselves.
    % 
    % It is tricky maintaining all this information but it does allow the
    % code to run much faster.

    % First extract the end nodes from each edgelist.  The nodes are numbered
    % so that the start node has number 2*edgenumber-1 and the end node has
    % number 2*edgenumber
    node = zeros(Nnodes, 2);
    for n = 1:Nedges
        node(2*n-1,:) = edgelist{n}(1,:);
        node(2*n  ,:) = edgelist{n}(end,:);     
    end
    
    % Now build the adjacency/connection matrices. 
    A = zeros(Nnodes);   % Adjacency matrix for nodes
    B = zeros(Nedges);   % Adjacency matrix for edges
    
    for n = 1:Nnodes-1
        for m = n+1:Nnodes
            % If nodes m & n are connected
            A(n,m) =  node(n,1)==node(m,1) && node(n,2)==node(m,2);
            A(m,n) = A(n,m);
            
            if A(n,m)
                edgen = fix((n+1)/2);
                edgem = fix((m+1)/2);           
                B(edgen, edgem) = 1;
                B(edgem, edgen) = 1;            
            end
        end
    end

    % If we sum the columns of the adjacency matrix we get the number of
    % other edgelists that are connected to an edgelist
    Nconnections = sum(A);   % Connection count array for nodes
    Econnections = sum(B);   % Connection count array for edges
    

    % Check every edge to see if any of its ends are connected to just one edge.
    % This should not happen, but occasionally does due to a problem in
    % EDGELINK.  Here we simply merge it with the edge it is connected to.
    % Ultimately I want to be able to remove this block of code.
    % I think there are also some cases that are (still) not properly handled
    % by CLEANEDGELIST and there may be a case for repeating this block of
    % code at the end for another final cleanup pass
    for n = 1:Nedges
        if ~B(n,n) && ~isempty(edgelist{n}) % if edge is not connected to itself
            [spurdegree, spurnode, startnode, sconns, endnode, econns] = connectioninfo(n);
            if sconns == 1
                node2merge = find(A(startnode,:));
                mergenodes(node2merge,startnode);
            end
            
            if ~isempty(edgelist{n})   % If we have not removed this edge in
                                       % the code above check the other end.
                if econns == 1
                    node2merge = find(A(endnode,:));
                    mergenodes(node2merge,endnode);
                end         
            end
        end
    end
    
    
    % Now check every edgelist, if the edgelength is below the minimum length
    % check if we should remove it.

    if minlength > 0
        
    for n = 1:Nedges
        
        [spurdegree, spurnode] = connectioninfo(n);
        
        if ~isempty(edgelist{n}) && edgelistlength(edgelist{n}) < minlength  

            % Remove unconnected lists, or lists that are only connected to
            % themselves. 
            if ~Econnections(n) || (Econnections(n)==1 && B(n,n) == 1)
                removeedge(n);
            
            % Process edges that are spurs coming from a 3-way junction.
            elseif spurdegree == 2
                %fprintf('%d is a spur\n',n)  %%debug

                linkingedges = find(B(n,:));
                
                if length(linkingedges) == 1 % We have a loop with a spur
                                             % coming from the join in the
                                             % loop
                   % Just remove the spur, leaving the loop intact.
                   removeedge(n);  
                   
                else   % Check the other edges coming from this point. If any
                       % are also spurs make sure we remove the shortest one
                   spurs = n;
                   len = edgelistlength(edgelist{n});
                   for i = 1:length(linkingedges)
                       spurdegree = connectioninfo(linkingedges(i));
                       if spurdegree
                           spurs = [spurs linkingedges(i)];  
                           len = [len edgelistlength(edgelist{linkingedges(i)})];  
                       end
                   end

                   linkingedges = [linkingedges n];
                   
                   [mn,i] = min(len);
                   edge2delete = spurs(i);
                   [spurdegree, spurnode] = connectioninfo(edge2delete);

                   nodes2merge = find(A(spurnode,:));
                   
                   if length(nodes2merge) ~= 2
                       error('attempt to merge other than 2 nodes');
                   end
                   
                   removeedge(edge2delete);             
                   mergenodes(nodes2merge(1),nodes2merge(2))               
                   
                end 

            % Look for spurs coming from 4-way junctions that are below the minimum length
            elseif spurdegree == 3
                removeedge(n);    % Just remove it, no subsequent merging needed.
            end
        end
    end
    
    % Final cleanup of any new isolated edges that might have been created by
    % removing spurs.  An edge is isolated if it has no connections to other
    % edges, or is only connected to itself (in a loop).
    
    for n = 1:Nedges
        if ~isempty(edgelist{n}) && edgelistlength(edgelist{n}) < minlength  
            if ~Econnections(n) || (Econnections(n)==1 && B(n,n) == 1)
                removeedge(n);          
            end
        end
    end
    
    end % if minlength > 0
    
    % Run through the edgelist and extract out the non-empty lists
    m = 0;
    for n = 1:Nedges
       if ~isempty(edgelist{n})
           m = m+1;
           nedgelist{m} = edgelist{n};
       end
    end
        
    
%---------------------------------------------------------------------- 
% Internal function to merge 2 edgelists together at the specified nodes and
% perform the necessary updates to the edge adjacency and node adjacency
% matrices and the connection count arrays

function mergenodes(n1,n2)
    
    edge1 = fix((n1+1)/2);   % Indices of the edges associated with the nodes
    edge2 = fix((n2+1)/2);    

    % Get indices of nodes at each end of the two edges
    s1 = 2*edge1-1; e1 = 2*edge1;
    s2 = 2*edge2-1; e2 = 2*edge2;    
    
    if edge1==edge2
        % We should not get here, but somehow we occasionally do
        % fprintf('Nodes %d %d\n',n1,n2)    %% debug
        % warning('Attempt to merge an edge with itself')
        return
    end
    
    if ~A(n1,n2)
        error('Attempt to merge nodes that are not connected');
    end
    
    if mod(n1,2)  % node n1 is the start of edge1
        flipedge1 = 1;  % edge1 will need to be reversed in order to join edge2
    else 
        flipedge1 = 0; 
    end
    
    if mod(n2,2)  % node n2 is the start of edge2       
        flipedge2 = 0;
    else
        flipedge2 = 1;
    end
    
    % Join edgelists together - with appropriate reordering depending on which
    % end is connected to which.  The result is stored in edge1
    
    if  ~flipedge1 && ~flipedge2 
        edgelist{edge1} = [edgelist{edge1}; edgelist{edge2}];

        A(e1,:) = A(e2,:);      A(:,e1) = A(:,e2);
        Nconnections(e1) = Nconnections(e2);
        
    elseif  ~flipedge1 && flipedge2
        edgelist{edge1} = [edgelist{edge1}; flipud(edgelist{edge2})];   
        
        A(e1,:) = A(s2,:);      A(:,e1) = A(:,s2);
        Nconnections(e1) = Nconnections(s2);

    elseif  flipedge1 && ~flipedge2
        edgelist{edge1} = [flipud(edgelist{edge1}); edgelist{edge2}]; 
        
        A(s1,:) = A(e1,:);      A(:,s1) = A(:,e1);
        A(e1,:) = A(e2,:);      A(:,e1) = A(:,e2);      
        Nconnections(s1) = Nconnections(e1);
        Nconnections(e1) = Nconnections(e2);    

    elseif  flipedge1 && flipedge2
        edgelist{edge1} = [flipud(edgelist{edge1}); flipud(edgelist{edge2})];

        A(s1,:) = A(e1,:);      A(:,s1) = A(:,e1);          
        A(e1,:) = A(s2,:);      A(:,e1) = A(:,s2);
        Nconnections(s1) = Nconnections(e1);
        Nconnections(e1) = Nconnections(s2);
        
    else
        fprintf('merging edges %d and %d\n',edge1, edge2); %%debug      
        error('We should not have got here - edgelists cannot be merged');
    end
    
    % Now correct the edge adjacency matrix to reflect the new arrangement
    % The edges that the new edge1 is connected to is all the edges that
    % edge1 and edge2 were connected to
    B(edge1,:) = B(edge1,:) | B(edge2,:);
    B(:,edge1) = B(:,edge1) | B(:,edge2);    
    B(edge1, edge1) = 0;

    % Recompute connection counts because we have shuffled the adjacency matrices
    Econnections = sum(B);
    Nconnections = sum(A);
  
    removeedge(edge2);  % Finally discard edge2
    
end  % end of mergenodes

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

% Function to provide information about the connections at each end of an
% edgelist 
%
%   [spurdegree, spurnode, startnode, sconns, endnode, econns] = connectioninfo(n)
%
%  spurdegree - If this is non-zero it indicates this edgelist is a spur, the
%               value is the number of edges this spur is connected to.
%  spurnode   - If this is a spur spurnode is the index of the node that is
%               connected to other edges, 0 otherwise.
%  startnode  - index of starting node of edgelist.
%  endnode    - index of end node of edgelist.
%  sconns     - number of connections to start node.
%  econns     - number of connections to end node.
    
function [spurdegree, spurnode, startnode, sconns, endnode, econns] = connectioninfo(n)

    if isempty(edgelist{n})
        spurdegree = 0; spurnode = 0;
        startnode = 0; sconns = 0; endnode = 0; econns = 0;
        return
    end
    
    startnode = 2*n-1;
    endnode   = 2*n;
    sconns = Nconnections(startnode);  % No of connections to start node
    econns = Nconnections(endnode);    % No of connections to end node    
    
    if sconns == 0 && econns >= 1
        spurdegree = econns;
        spurnode = endnode;
    elseif sconns >= 1 && econns == 0
        spurdegree = sconns;
        spurnode = startnode;   
    else
        spurdegree = 0;
        spurnode = 0;
    end
    
end

%--------------------------------------------------------------------
% Function to remove an edgelist and perform the necessary updates to the edge
% adjacency and node adjacency matrices and the connection count arrays

function removeedge(n)
    
    edgelist{n} = [];
    Econnections = Econnections - B(n,:); 
    Econnections(n) = 0; 
    B(n,:) = 0;
    B(:,n) = 0;
    
    nodes2delete = [2*n-1, 2*n];
    
    Nconnections = Nconnections - A(nodes2delete(1),:);    
    Nconnections = Nconnections - A(nodes2delete(2),:);        
    
    A(nodes2delete, :) = 0;
    A(:, nodes2delete) = 0;                 
    
end

%--------------------------------------------------------------------
% Function to compute the path length of an edgelist

function l = edgelistlength(edgelist)
    l = sum(sqrt(sum((edgelist(1:end-1,:)-edgelist(2:end,:)).^2, 2)));
end

%--------------------------------------------------------------------
end % End of cleanedgelists