intersectLines

INTERSECTLINES Return all intersection points of N lines in 2D.

   PT = intersectLines(L1, L2);
   returns the intersection point of lines L1 and L2. L1 and L2 are 1-by-4
   row arrays, containing parametric representation of each line (in the
   form [x0 y0 dx dy], see 'createLine' for details).
   
   In case of colinear lines, returns [Inf Inf].
   In case of parallel but not colinear lines, returns [NaN NaN].

   If each input is [N*4] array, the result is a [N*2] array containing
   intersections of each couple of lines.
   If one of the input has N rows and the other 1 row, the result is a
   [N*2] array.

   PT = intersectLines(L1, L2, EPS);
   Specifies the tolerance for detecting parallel lines. Default is 1e-14.

   Example
   line1 = createLine([0 0], [10 10]);
   line2 = createLine([0 10], [10 0]);
   point = intersectLines(line1, line2)
   point = 
       5   5

   See also
   lines2d, edges2d, intersectEdges, intersectLineEdge
   intersectLineCircle

0001 function point = intersectLines(line1, line2, varargin)
0002 %INTERSECTLINES Return all intersection points of N lines in 2D.
0003 %
0004 %   PT = intersectLines(L1, L2);
0005 %   returns the intersection point of lines L1 and L2. L1 and L2 are 1-by-4
0006 %   row arrays, containing parametric representation of each line (in the
0007 %   form [x0 y0 dx dy], see 'createLine' for details).
0008 %
0009 %   In case of colinear lines, returns [Inf Inf].
0010 %   In case of parallel but not colinear lines, returns [NaN NaN].
0011 %
0012 %   If each input is [N*4] array, the result is a [N*2] array containing
0013 %   intersections of each couple of lines.
0014 %   If one of the input has N rows and the other 1 row, the result is a
0015 %   [N*2] array.
0016 %
0017 %   PT = intersectLines(L1, L2, EPS);
0018 %   Specifies the tolerance for detecting parallel lines. Default is 1e-14.
0019 %
0020 %   Example
0021 %   line1 = createLine([0 0], [10 10]);
0022 %   line2 = createLine([0 10], [10 0]);
0023 %   point = intersectLines(line1, line2)
0024 %   point =
0025 %       5   5
0026 %
0027 %   See also
0028 %   lines2d, edges2d, intersectEdges, intersectLineEdge
0029 %   intersectLineCircle
0030 
0031 % ------
0032 % Author: David Legland
0033 % E-mail: david.legland@inrae.fr
0034 % Created: 2003-10-31
0035 % Copyright 2003-2024 INRA - TPV URPOI - BIA IMASTE
0036 
0037 %% Process input arguments
0038 
0039 % extract tolerance
0040 tol = 1e-14;
0041 if ~isempty(varargin)
0042     tol = varargin{1};
0043 end
0044 
0045 % check size of each input
0046 N1 = size(line1, 1);
0047 N2 = size(line2, 1);
0048 N = max(N1, N2);
0049 if N1 ~= N2 && N1*N2 ~= N
0050     error('matGeom:IntersectLines:IllegalArgument', ...
0051         'The two input arguments must have same number of lines');
0052 end
0053 
0054 
0055 %% Check parallel and colinear lines
0056 
0057 % coordinate differences of origin points
0058 dx = bsxfun(@minus, line2(:,1), line1(:,1));
0059 dy = bsxfun(@minus, line2(:,2), line1(:,2));
0060 
0061 % indices of parallel lines
0062 denom = line1(:,3) .* line2(:,4) - line2(:,3) .* line1(:,4);
0063 par = abs(denom) < tol;
0064 
0065 % indices of colinear lines
0066 col = abs(dx .* line1(:,4) - dy .* line1(:,3)) < tol & par ;
0067 
0068 % initialize result array
0069 x0 = zeros(N, 1);
0070 y0 = zeros(N, 1);
0071 
0072 % initialize result for parallel lines
0073 x0(col) = Inf;
0074 y0(col) = Inf;
0075 x0(par & ~col) = NaN;
0076 y0(par & ~col) = NaN;
0077 
0078 % in case all line couples are parallel, return
0079 if all(par)
0080     point = [x0 y0];
0081     return;
0082 end
0083 
0084 
0085 %% Extract coordinates of itnersecting lines
0086 
0087 % indices of intersecting lines
0088 inds = ~par;
0089 
0090 % extract base coordinates of first lines
0091 if N1 > 1
0092     line1 = line1(inds,:);
0093 end
0094 x1 =  line1(:,1);
0095 y1 =  line1(:,2);
0096 dx1 = line1(:,3);
0097 dy1 = line1(:,4);
0098 
0099 % extract base coordinates of second lines
0100 if N2 > 1
0101     line2 = line2(inds,:);
0102 end
0103 x2 =  line2(:,1);
0104 y2 =  line2(:,2);
0105 dx2 = line2(:,3);
0106 dy2 = line2(:,4);
0107 
0108 % re-compute coordinate differences of origin points
0109 dx = bsxfun(@minus, line2(:,1), line1(:,1));
0110 dy = bsxfun(@minus, line2(:,2), line1(:,2));
0111 
0112 
0113 %% Compute intersection points
0114 
0115 denom = denom(inds);
0116 x0(inds) = (x2 .* dy2 .* dx1 - dy .* dx1 .* dx2 - x1 .* dy1 .* dx2) ./ denom ;
0117 y0(inds) = (dx .* dy1 .* dy2 + y1 .* dx1 .* dy2 - y2 .* dx2 .* dy1) ./ denom ;
0118 
0119 % concatenate result
0120 point = [x0 y0];