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

   ---------
   author : David Legland
   INRA - TPV URPOI - BIA IMASTE
   created the 31/10/2003.

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 %   INRA - TPV URPOI - BIA IMASTE
0034 %   created the 31/10/2003.
0035 %
0036 
0037 %   HISTORY
0038 %   2004-02-19 add support for multiple lines.
0039 %   2007-03-08 update doc
0040 %   2011-10-07 code cleanup
0041 
0042 
0043 %% Process input arguments
0044 
0045 % extract tolerance
0046 tol = 1e-14;
0047 if ~isempty(varargin)
0048     tol = varargin{1};
0049 end
0050 
0051 % check size of each input
0052 N1 = size(line1, 1);
0053 N2 = size(line2, 1);
0054 N = max(N1, N2);
0055 if N1 ~= N2 && N1*N2 ~= N
0056     error('matGeom:IntersectLines:IllegalArgument', ...
0057         'The two input arguments must have same number of lines');
0058 end
0059 
0060 
0061 %% Check parallel and colinear lines
0062 
0063 % coordinate differences of origin points
0064 dx = bsxfun(@minus, line2(:,1), line1(:,1));
0065 dy = bsxfun(@minus, line2(:,2), line1(:,2));
0066 
0067 % indices of parallel lines
0068 denom = line1(:,3) .* line2(:,4) - line2(:,3) .* line1(:,4);
0069 par = abs(denom) < tol;
0070 
0071 % indices of colinear lines
0072 col = abs(dx .* line1(:,4) - dy .* line1(:,3)) < tol & par ;
0073 
0074 % initialize result array
0075 x0 = zeros(N, 1);
0076 y0 = zeros(N, 1);
0077 
0078 % initialize result for parallel lines
0079 x0(col) = Inf;
0080 y0(col) = Inf;
0081 x0(par & ~col) = NaN;
0082 y0(par & ~col) = NaN;
0083 
0084 % in case all line couples are parallel, return
0085 if all(par)
0086     point = [x0 y0];
0087     return;
0088 end
0089 
0090 
0091 %% Extract coordinates of itnersecting lines
0092 
0093 % indices of intersecting lines
0094 inds = ~par;
0095 
0096 % extract base coordinates of first lines
0097 if N1 > 1
0098     line1 = line1(inds,:);
0099 end
0100 x1 =  line1(:,1);
0101 y1 =  line1(:,2);
0102 dx1 = line1(:,3);
0103 dy1 = line1(:,4);
0104 
0105 % extract base coordinates of second lines
0106 if N2 > 1
0107     line2 = line2(inds,:);
0108 end
0109 x2 =  line2(:,1);
0110 y2 =  line2(:,2);
0111 dx2 = line2(:,3);
0112 dy2 = line2(:,4);
0113 
0114 % re-compute coordinate differences of origin points
0115 dx = bsxfun(@minus, line2(:,1), line1(:,1));
0116 dy = bsxfun(@minus, line2(:,2), line1(:,2));
0117 
0118 
0119 %% Compute intersection points
0120 
0121 denom = denom(inds);
0122 x0(inds) = (x2 .* dy2 .* dx1 - dy .* dx1 .* dx2 - x1 .* dy1 .* dx2) ./ denom ;
0123 y0(inds) = (dx .* dy1 .* dy2 + y1 .* dx1 .* dy2 - y2 .* dx2 .* dy1) ./ denom ;
0124 
0125 % concatenate result
0126 point = [x0 y0];