function new_edges = splitEdges( edges )
% splitEdges: split line segments (=edges) into smaller ones based on 
%             intersections between them
% 
% isIntersect.m, intersectEdges.m, isPointOnEdge.m are used.
% isIntersect.m - algorithm (www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/)
% intersectEdges.m, isPointOnEdge.m are written by David Legland (https://github.com/mattools/matGeom)
%
% input:
%	edges - N-by-4 array, storing N line segments. 
%    		Each line segment (edges(i,:)) is given by coordinates 
%           of two points : [x1 y1 x2 y2]. 
%
% example:
%
%   edges = [ 0 0 5 0; 2 0 4 0; 1 5 1 -3; 1 5 3 -2 ];
%   new_edges = splitEdges( edges );
%   drawEdge( new_edges, 'linewidth', 1.5 );
%
% Written by Jiexian Ma at SUNY Binghamton (mjx0799@gmail.com)
%
       
    % find intersection between edges
    pi_cell = edgexedge( edges );
    % sort intersections (pi_cell) on individual edge
    pi_cell = sortPoint( pi_cell );
    
    % sort 2 node in each edge
    edges = sortNode( edges );
    % break up each edge by intersections on them
    new_edges = breakEdges( edges, pi_cell );

end


function pi_cell = edgexedge( edges )
% find intersection between edges
% pi_cell - cell column vector, record intersections
% pi_cell{i} - n-by-2 array, save x y data of intersections on edges(i,:)

    % point of intersection between edges
    pi_cell = cell( size(edges,1), 1 );  % cell array
    % pi_cell{i} - n by 2 array, save x y data

    TOLERANCE = 10^(-10);   % for ismembertol

    % find intersection between edges
    for i = 1: size(edges,1)-1
        for j = i+1: size(edges,1)

            e1 = edges( i, : );
            e2 = edges( j, : );
            
            if ~ isIntersect( e1, e2 )
                % case of no intersect by isIntersect.m
                xyi = [NaN NaN];
            else
                % find intersection between two edges using
                % intersectEdges.m
                % 0-2 points
                xyi = intersectEdges( e1, e2, TOLERANCE );
            end
                
            node1 = e1( [1 2] );
            node2 = e1( [3 4] );
            node3 = e2( [1 2] );
            node4 = e2( [3 4] );
            
            if all( isnan(xyi) )	% case of no intersect: xyi = [NaN NaN]
                % do nothing
            else                    % case of intersect
                
                if all( isinf(xyi) ) % case of intersect but colinear edges
                                     % xyi = [inf inf]
                    % find out intersects (colinear case)
                    % using isPointOnEdge.m
                    xytemp = [];
                    if isPointOnEdge( node1, e2, TOLERANCE )
                        xytemp = [ xytemp; node1];
                    end
                    if isPointOnEdge( node2, e2, TOLERANCE )
                        xytemp = [ xytemp; node2];
                    end
                    if isPointOnEdge( node3, e1, TOLERANCE )
                        xytemp = [ xytemp; node3];
                    end
                    if isPointOnEdge( node4, e1, TOLERANCE )
                        xytemp = [ xytemp; node4];
                    end
                    xyi = xytemp;
                end
                
                % exclude the original nodes, and save xyi
                % compare xyi with the nodes of edge1
                ind1 = ~ismembertol( xyi, [node1;node2], TOLERANCE,'ByRows',true);
                % store non-node xyi
                pi_cell{i} = [ pi_cell{i}; xyi( ind1, : ) ];

                % compare xyi with the nodes of edge2
                ind2 = ~ismembertol( xyi, [node3;node4], TOLERANCE,'ByRows',true);
                % store non-node xyi
                pi_cell{j} = [ pi_cell{j}; xyi( ind2, : ) ];
            end

        end
    end
end

function pi_cell = sortPoint( pi_cell )
% sort intersections in pi_cell

    for i = 1: length( pi_cell )
        temp = unique( pi_cell{i}, 'rows' );
        
        % sort intersections
        if size( temp, 1 ) > 1
            if temp(1,1) ~= temp(2,1)
                % sort points according to x
                [ ~, ind ] = sort( temp(:,1) );
                temp = temp( ind, : );
            else
                % sort points according to y
                [ ~, ind ] = sort( temp(:,2) );
                temp = temp( ind, : );
            end
        end
        
        pi_cell{i} = temp;
    end
end

function edges = sortNode( edges )
% sort 2 node in each edge

    for i = 1: size( edges, 1 )
        temp = [ edges(i,1:2); edges(i,3:4)];

        % sort
        if temp(1,1) ~= temp(2,1)
            % sort points according to x
            [ ~, ind ] = sort( temp(:,1) );
            temp = temp( ind, : );
        else
            % sort points according to y
            [ ~, ind ] = sort( temp(:,2) );
            temp = temp( ind, : );
        end

        edges( i, : ) = [ temp(1,:) temp(2,:) ];
    end
end

function new_edges = breakEdges( edges, pi_cell )
% break up edges

    % total number of edges after breaking up
    sum_len = 0;
    for i = 1: length( pi_cell )
        sum_len = sum_len + size( pi_cell{i}, 1 ) + 1;
    end

    new_edges = zeros( sum_len, 4 );
    k = 1;  % index for new_edges

    for i = 1: length( pi_cell )
        if size( pi_cell{i}, 1 ) == 0
            % no intection on this edge
            new_edges(k,:) = edges(i,:);  % do nothing
            k = k + 1;
            
        elseif size( pi_cell{i}, 1 ) == 1
            % 1 intection on this edge
            new_edges(k,:) = [ edges(i,1:2) pi_cell{i}(1,1:2) ];
            k = k + 1;

            new_edges(k,:) = [ pi_cell{i}(end,1:2) edges(i,3:4) ];
            k = k + 1;

        elseif size( pi_cell{i}, 1 ) > 1
            % >1 intection on this edge
            for j = 1: size( pi_cell{i}, 1 ) - 1
                new_edges(k,:) = [ pi_cell{i}(j,1:2) pi_cell{i}(j+1,1:2) ];
                k = k + 1;
            end

            new_edges(k,:) = [ edges(i,1:2) pi_cell{i}(1,1:2) ];
            k = k + 1;

            new_edges(k,:) = [ pi_cell{i}(end,1:2) edges(i,3:4) ];
            k = k + 1;
        end
    end
end