lineSegmentIntersect

LINESEGMENTINTERSECT Intersections of line segments.
   OUT = LINESEGMENTINTERSECT(XY1,XY2) finds the 2D Cartesian Coordinates of
   intersection points between the set of line segments given in XY1 and XY2.

   XY1 and XY2 are N1x4 and N2x4 matrices. Rows correspond to line segments. 
   Each row is of the form [x1 y1 x2 y2] where (x1,y1) is the start point and 
   (x2,y2) is the end point of a line segment:

                  Line Segment
       o--------------------------------o
       ^                                ^
    (x1,y1)                          (x2,y2)

   OUT is a structure with fields:

   'intAdjacencyMatrix' : N1xN2 indicator matrix where the entry (i,j) is 1 if
       line segments XY1(i,:) and XY2(j,:) intersect.

   'intMatrixX' : N1xN2 matrix where the entry (i,j) is the X coordinate of the
       intersection point between line segments XY1(i,:) and XY2(j,:).

   'intMatrixY' : N1xN2 matrix where the entry (i,j) is the Y coordinate of the
       intersection point between line segments XY1(i,:) and XY2(j,:).

   'intNormalizedDistance1To2' : N1xN2 matrix where the (i,j) entry is the
       normalized distance from the start point of line segment XY1(i,:) to the
       intersection point with XY2(j,:).

   'intNormalizedDistance2To1' : N1xN2 matrix where the (i,j) entry is the
       normalized distance from the start point of line segment XY1(j,:) to the
       intersection point with XY2(i,:).

   'parAdjacencyMatrix' : N1xN2 indicator matrix where the (i,j) entry is 1 if
       line segments XY1(i,:) and XY2(j,:) are parallel.

   'coincAdjacencyMatrix' : N1xN2 indicator matrix where the (i,j) entry is 1 
       if line segments XY1(i,:) and XY2(j,:) are coincident.

0001 function out = lineSegmentIntersect(XY1,XY2)
0002 %LINESEGMENTINTERSECT Intersections of line segments.
0003 %   OUT = LINESEGMENTINTERSECT(XY1,XY2) finds the 2D Cartesian Coordinates of
0004 %   intersection points between the set of line segments given in XY1 and XY2.
0005 %
0006 %   XY1 and XY2 are N1x4 and N2x4 matrices. Rows correspond to line segments.
0007 %   Each row is of the form [x1 y1 x2 y2] where (x1,y1) is the start point and
0008 %   (x2,y2) is the end point of a line segment:
0009 %
0010 %                  Line Segment
0011 %       o--------------------------------o
0012 %       ^                                ^
0013 %    (x1,y1)                          (x2,y2)
0014 %
0015 %   OUT is a structure with fields:
0016 %
0017 %   'intAdjacencyMatrix' : N1xN2 indicator matrix where the entry (i,j) is 1 if
0018 %       line segments XY1(i,:) and XY2(j,:) intersect.
0019 %
0020 %   'intMatrixX' : N1xN2 matrix where the entry (i,j) is the X coordinate of the
0021 %       intersection point between line segments XY1(i,:) and XY2(j,:).
0022 %
0023 %   'intMatrixY' : N1xN2 matrix where the entry (i,j) is the Y coordinate of the
0024 %       intersection point between line segments XY1(i,:) and XY2(j,:).
0025 %
0026 %   'intNormalizedDistance1To2' : N1xN2 matrix where the (i,j) entry is the
0027 %       normalized distance from the start point of line segment XY1(i,:) to the
0028 %       intersection point with XY2(j,:).
0029 %
0030 %   'intNormalizedDistance2To1' : N1xN2 matrix where the (i,j) entry is the
0031 %       normalized distance from the start point of line segment XY1(j,:) to the
0032 %       intersection point with XY2(i,:).
0033 %
0034 %   'parAdjacencyMatrix' : N1xN2 indicator matrix where the (i,j) entry is 1 if
0035 %       line segments XY1(i,:) and XY2(j,:) are parallel.
0036 %
0037 %   'coincAdjacencyMatrix' : N1xN2 indicator matrix where the (i,j) entry is 1
0038 %       if line segments XY1(i,:) and XY2(j,:) are coincident.
0039 
0040 % Version: 1.00, April 03, 2010
0041 % Version: 1.10, April 10, 2010
0042 % Author:  U. Murat Erdem
0043 
0044 % CHANGELOG:
0045 %
0046 % Ver. 1.00:
0047 %   -Initial release.
0048 %
0049 % Ver. 1.10:
0050 %   - Changed the input parameters. Now the function accepts two sets of line
0051 %   segments. The intersection analysis is done between these sets and not in
0052 %   the same set.
0053 %   - Changed and added fields of the output. Now the analysis provides more
0054 %   information about the intersections and line segments.
0055 %   - Performance tweaks.
0056 
0057 % I opted not to call this 'curve intersect' because it would be misleading
0058 % unless you accept that curves are pairwise linear constructs.
0059 % I tried to put emphasis on speed by vectorizing the code as much as possible.
0060 % There should still be enough room to optimize the code but I left those out
0061 % for the sake of clarity.
0062 % The math behind is given in:
0063 %   http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline2d/
0064 % If you really are interested in squeezing as much horse power as possible out
0065 % of this code I would advise to remove the argument checks and tweak the
0066 % creation of the OUT a little bit.
0067 
0068 %%% Argument check.
0069 %-------------------------------------------------------------------------------
0070 
0071 validateattributes(XY1,{'numeric'},{'2d','finite'});
0072 validateattributes(XY2,{'numeric'},{'2d','finite'});
0073 
0074 [n_rows_1,n_cols_1] = size(XY1);
0075 [n_rows_2,n_cols_2] = size(XY2);
0076 
0077 if n_cols_1 ~= 4 || n_cols_2 ~= 4
0078     error('Arguments must be a Nx4 matrices.');
0079 end
0080 
0081 %%% Prepare matrices for vectorized computation of line intersection points.
0082 %-------------------------------------------------------------------------------
0083 X1 = repmat(XY1(:,1),1,n_rows_2);
0084 X2 = repmat(XY1(:,3),1,n_rows_2);
0085 Y1 = repmat(XY1(:,2),1,n_rows_2);
0086 Y2 = repmat(XY1(:,4),1,n_rows_2);
0087 
0088 XY2 = XY2';
0089 
0090 X3 = repmat(XY2(1,:),n_rows_1,1);
0091 X4 = repmat(XY2(3,:),n_rows_1,1);
0092 Y3 = repmat(XY2(2,:),n_rows_1,1);
0093 Y4 = repmat(XY2(4,:),n_rows_1,1);
0094 
0095 X4_X3 = (X4-X3);
0096 Y1_Y3 = (Y1-Y3);
0097 Y4_Y3 = (Y4-Y3);
0098 X1_X3 = (X1-X3);
0099 X2_X1 = (X2-X1);
0100 Y2_Y1 = (Y2-Y1);
0101 
0102 numerator_a = X4_X3 .* Y1_Y3 - Y4_Y3 .* X1_X3;
0103 numerator_b = X2_X1 .* Y1_Y3 - Y2_Y1 .* X1_X3;
0104 denominator = Y4_Y3 .* X2_X1 - X4_X3 .* Y2_Y1;
0105 
0106 u_a = numerator_a ./ denominator;
0107 u_b = numerator_b ./ denominator;
0108 
0109 % Find the adjacency matrix A of intersecting lines.
0110 INT_X = X1+X2_X1.*u_a;
0111 INT_Y = Y1+Y2_Y1.*u_a;
0112 INT_B = (u_a >= 0) & (u_a <= 1) & (u_b >= 0) & (u_b <= 1);
0113 PAR_B = denominator == 0;
0114 COINC_B = (numerator_a == 0 & numerator_b == 0 & PAR_B);
0115 
0116 
0117 % Arrange output.
0118 out.intAdjacencyMatrix = INT_B;
0119 out.intMatrixX = INT_X .* INT_B;
0120 out.intMatrixY = INT_Y .* INT_B;
0121 out.intNormalizedDistance1To2 = u_a;
0122 out.intNormalizedDistance2To1 = u_b;
0123 out.parAdjacencyMatrix = PAR_B;
0124 out.coincAdjacencyMatrix= COINC_B;
0125 
0126 end
0127