intersectPolylines

INTERSECTPOLYLINES Find the common points between 2 polylines.

   INTERS = intersectPolylines(POLY1, POLY2)
   Returns the intersection points between two polylines. Each polyline is
   defined by a N-by-2 array representing coordinates of its vertices: 
   [X1 Y1 ; X2 Y2 ; ... ; XN YN]
   INTERS is a NP-by-2 array containing coordinates of intersection
   points.

   INTERS = intersectPolylines(POLY1)
   Compute self-intersections of the polyline.

   Example
   % Compute intersection points between 2 simple polylines
     poly1 = [20 10 ; 20 50 ; 60 50 ; 60 10];
     poly2 = [10 40 ; 30 40 ; 30 60 ; 50 60 ; 50 40 ; 70 40];
     pts = intersectPolylines(poly1, poly2);
     figure; hold on; 
     drawPolyline(poly1, 'b');
     drawPolyline(poly2, 'm');
     drawPoint(pts);
     axis([0 80 0 80]);

   This function is largely based on the 'interX' function, found on the
   FileExchange:
   https://fr.mathworks.com/matlabcentral/fileexchange/22441-curve-intersections
   
   See also
   polygons2d, polylineSelfIntersections, intersectLinePolygon

0001 function pts = intersectPolylines(poly1, varargin)
0002 %INTERSECTPOLYLINES Find the common points between 2 polylines.
0003 %
0004 %   INTERS = intersectPolylines(POLY1, POLY2)
0005 %   Returns the intersection points between two polylines. Each polyline is
0006 %   defined by a N-by-2 array representing coordinates of its vertices:
0007 %   [X1 Y1 ; X2 Y2 ; ... ; XN YN]
0008 %   INTERS is a NP-by-2 array containing coordinates of intersection
0009 %   points.
0010 %
0011 %   INTERS = intersectPolylines(POLY1)
0012 %   Compute self-intersections of the polyline.
0013 %
0014 %   Example
0015 %   % Compute intersection points between 2 simple polylines
0016 %     poly1 = [20 10 ; 20 50 ; 60 50 ; 60 10];
0017 %     poly2 = [10 40 ; 30 40 ; 30 60 ; 50 60 ; 50 40 ; 70 40];
0018 %     pts = intersectPolylines(poly1, poly2);
0019 %     figure; hold on;
0020 %     drawPolyline(poly1, 'b');
0021 %     drawPolyline(poly2, 'm');
0022 %     drawPoint(pts);
0023 %     axis([0 80 0 80]);
0024 %
0025 %   This function is largely based on the 'interX' function, found on the
0026 %   FileExchange:
0027 %   https://fr.mathworks.com/matlabcentral/fileexchange/22441-curve-intersections
0028 %
0029 %   See also
0030 %   polygons2d, polylineSelfIntersections, intersectLinePolygon
0031 %
0032 
0033 % ------
0034 % Author: David Legland
0035 % e-mail: david.legland@inra.fr
0036 % Created: 2009-06-15,    using Matlab 7.7.0.471 (R2008b)
0037 % Copyright 2009 INRA - Cepia Software Platform.
0038 
0039 % The code is a slight rewritting of the interX function, consisting in
0040 % avoiding argument transposition in the begining of the function. Comment
0041 % of original submission are kept here for information.
0042 %
0043 %INTERX Intersection of curves
0044 %   P = INTERX(L1,L2) returns the intersection points of two curves L1
0045 %   and L2. The curves L1,L2 can be either closed or open and are described
0046 %   by two-row-matrices, where each row contains its x- and y- coordinates.
0047 %   The intersection of groups of curves (e.g. contour lines, multiply
0048 %   connected regions etc) can also be computed by separating them with a
0049 %   column of NaNs as for example
0050 %
0051 %         L  = [x11 x12 x13 ... NaN x21 x22 x23 ...;
0052 %               y11 y12 y13 ... NaN y21 y22 y23 ...]
0053 %
0054 %   P has the same structure as L1 and L2, and its rows correspond to the
0055 %   x- and y- coordinates of the intersection points of L1 and L2. If no
0056 %   intersections are found, the returned P is empty.
0057 %
0058 %   P = INTERX(L1) returns the self-intersection points of L1. To keep
0059 %   the code simple, the points at which the curve is tangent to itself are
0060 %   not included. P = INTERX(L1,L1) returns all the points of the curve
0061 %   together with any self-intersection points.
0062 %
0063 %   Example:
0064 %       t = linspace(0,2*pi);
0065 %       r1 = sin(4*t)+2;  x1 = r1.*cos(t); y1 = r1.*sin(t);
0066 %       r2 = sin(8*t)+2;  x2 = r2.*cos(t); y2 = r2.*sin(t);
0067 %       P = InterX([x1;y1],[x2;y2]);
0068 %       plot(x1,y1,x2,y2,P(1,:),P(2,:),'ro')
0069 %
0070 %   Author : NS
0071 %   Version: 3.0, 21 Sept. 2010
0072 %
0073 %   Two words about the algorithm: Most of the code is self-explanatory.
0074 %   The only trick lies in the calculation of C1 and C2. To be brief, this
0075 %   is essentially the two-dimensional analog of the condition that needs
0076 %   to be satisfied by a function F(x) that has a zero in the interval
0077 %   [a,b], namely
0078 %           F(a)*F(b) <= 0
0079 %   C1 and C2 exactly do this for each segment of curves 1 and 2
0080 %   respectively. If this condition is satisfied simultaneously for two
0081 %   segments then we know that they will cross at some point.
0082 %   Each factor of the 'C' arrays is essentially a matrix containing
0083 %   the numerators of the signed distances between points of one curve
0084 %   and line segments of the other.
0085 
0086 
0087 % Check number of inputs
0088 narginchk(1, 2);
0089 
0090 % Specific init depending on number of inputs
0091 if nargin == 1
0092     % Compute self-intersections
0093     % -> Avoid the inclusion of common points
0094     poly2 = poly1;
0095     hF = @lt;
0096 else
0097     % Compute intersections between distinct lines
0098     poly2 = varargin{1}; 
0099     hF = @le;
0100 end
0101 
0102 % Get coordinates of polyline vertices
0103 x1 = poly1(:,1);  
0104 x2 = poly2(:,1)';
0105 y1 = poly1(:,2);  
0106 y2 = poly2(:,2)';
0107 
0108 % differentiate coordinate arrays
0109 dx1 = diff(x1); dy1 = diff(y1);
0110 dx2 = diff(x2); dy2 = diff(y2);
0111 
0112 % Determine 'signed distances'
0113 S1 = dx1 .* y1(1:end-1) - dy1 .* x1(1:end-1);
0114 S2 = dx2 .* y2(1:end-1) - dy2 .* x2(1:end-1);
0115 
0116 C1 = feval(hF, D(bsxfun(@times,dx1,y2) - bsxfun(@times,dy1,x2), S1), 0);
0117 C2 = feval(hF, D((bsxfun(@times,y1,dx2) - bsxfun(@times,x1,dy2))', S2'), 0)';
0118 
0119 % Obtain the segments where an intersection is expected
0120 [i, j] = find(C1 & C2); 
0121 
0122 % Process case of no intersection
0123 if isempty(i)
0124     pts = zeros(0, 2);
0125     return;
0126 end
0127 
0128 % Transpose and prepare for output
0129 i=i'; dx2=dx2'; dy2=dy2'; S2 = S2';
0130 L = dy2(j).*dx1(i) - dy1(i).*dx2(j);
0131 
0132 % Avoid divisions by zero
0133 i = i(L~=0);
0134 j = j(L~=0);
0135 L = L(L~=0);
0136 
0137 % Solve system of eqs to get the common points
0138 res = [dx2(j).*S1(i) - dx1(i).*S2(j), dy2(j).*S1(i) - dy1(i).*S2(j)] ./ [L L];
0139 pts = unique(res, 'rows');
0140 
0141 % Innre function computing a kind of cross-product
0142 function u = D(x,y)
0143     u = bsxfun(@minus, x(:,1:end-1), y) .* bsxfun(@minus, x(:,2:end), y);
0144 end
0145 
0146 end