


INTSECL - Intersection coordinates of two line segments.
[XI,YI] = INTSECL(X1,Y1,X2,Y2) where all 4
arguments are 2 by N matrices with coordinates
of ends of N pairs of line segments (so that
the command such as PLOT(X1,Y1,X2,Y2) will plot
these pairs of lines).
Returns 1 by N vectors XI and YI consisting of
coordinates of intersection points of each of N
pairs of lines.
Special cases:
When a line segment is degenerate into a point
and does not lie on line through the other segment
of a pair returns XI=NaN while YI has the following
values: 1 - when the first segment in a pair is
degenerate, 2 - second segment, 0 - both segments
are degenerate.
When a pair of line segments is parallel, returns
XI = Inf while YI is 1 for coincident lines,
0 - for parallel non-coincident ones.
INTSECL(X1,Y1,X2,Y2,TOL) also specifies tolerance
in detecting coincident points in different line
segments.
Copyright (c) 1995 by Kirill K. Pankratov
kirill@plume.mit.edu
04/15/94, 08/14/94, 05/10/95, 08/23/95

0001 function [xo,yo] = intsecl(x1,y1,x2,y2,tol) 0002 %INTSECL - Intersection coordinates of two line segments. 0003 % 0004 % [XI,YI] = INTSECL(X1,Y1,X2,Y2) where all 4 0005 % arguments are 2 by N matrices with coordinates 0006 % of ends of N pairs of line segments (so that 0007 % the command such as PLOT(X1,Y1,X2,Y2) will plot 0008 % these pairs of lines). 0009 % Returns 1 by N vectors XI and YI consisting of 0010 % coordinates of intersection points of each of N 0011 % pairs of lines. 0012 % 0013 % Special cases: 0014 % When a line segment is degenerate into a point 0015 % and does not lie on line through the other segment 0016 % of a pair returns XI=NaN while YI has the following 0017 % values: 1 - when the first segment in a pair is 0018 % degenerate, 2 - second segment, 0 - both segments 0019 % are degenerate. 0020 % When a pair of line segments is parallel, returns 0021 % XI = Inf while YI is 1 for coincident lines, 0022 % 0 - for parallel non-coincident ones. 0023 % INTSECL(X1,Y1,X2,Y2,TOL) also specifies tolerance 0024 % in detecting coincident points in different line 0025 % segments. 0026 % 0027 % Copyright (c) 1995 by Kirill K. Pankratov 0028 % kirill@plume.mit.edu 0029 % 04/15/94, 08/14/94, 05/10/95, 08/23/95 0030 0031 0032 % Defaults and parameters ......................... 0033 tol_dflt = 0; % Tolerance for coincident points 0034 is_chk = 1; % Check input arguments 0035 0036 % Handle input .................................... 0037 if nargin==0, help intsecl, return, end 0038 if nargin<4 % Check if all 4 entered 0039 error(' Not enough input arguments') 0040 end 0041 if nargin<5, tol = tol_dflt; end 0042 if tol < 0, is_chk = 0; tol = 0; end 0043 0044 % Check the format of arguments ....... 0045 if is_chk 0046 [x1,y1,x2,y2] = linechk(x1,y1,x2,y2); 0047 end 0048 0049 0050 % Auxillary 0051 o2 = ones(2,1); 0052 i_pt1 = []; i_pt2 = []; i_pt12 = []; 0053 0054 % Make first points origins ........... 0055 xo = x1(1,:); 0056 yo = y1(1,:); 0057 x2 = x2-xo(o2,:); 0058 y2 = y2-yo(o2,:); 0059 0060 % Differences of first segments ....... 0061 a = x1(2,:)-x1(1,:); 0062 b = y1(2,:)-y1(1,:); 0063 s = sqrt(a.^2+b.^2); % Lengths of first segments 0064 i_pt1 = find(~s); 0065 s(i_pt1) = ones(size(i_pt1)); 0066 rr = rand(size(i_pt1)); 0067 a(i_pt1) = cos(rr); 0068 b(i_pt1) = sin(rr); 0069 0070 % Normalize by length ................. 0071 a = a./s; b = b./s; 0072 0073 % Rotate coordinates of the second segment 0074 tmp = x2.*a(o2,:)+y2.*b(o2,:); 0075 y2 = -x2.*b(o2,:)+y2.*a(o2,:); 0076 x2 = tmp; 0077 0078 % Calculate differences in second segments 0079 s = x2(2,:)-x2(1,:); 0080 tmp = y2(2,:)-y2(1,:); 0081 cc = tmp(i_pt1); 0082 0083 % Find some degenerate cases ....................... 0084 0085 % Find zeros in differences 0086 izy2 = find(~tmp); 0087 tmp(izy2) = ones(size(izy2)); 0088 0089 % Find degenerate and parallel segments 0090 bool = ~s(izy2); 0091 i_par = izy2(~bool); 0092 i_pt2 = izy2(bool); 0093 0094 bool = abs(y2(1,i_pt2))<=tol; 0095 i_pt2_off = i_pt2(~bool); 0096 i_pt2_on = i_pt2(bool); 0097 0098 if ~isempty(i_par) 0099 bool = abs(y2(1,i_par))<=tol; 0100 i_par_off = i_par(~bool); 0101 i_par_on = i_par(bool); 0102 end 0103 0104 % Calculate intercept with rotated x-axis .......... 0105 tmp = s./tmp; % Slope 0106 tmp = x2(1,:)-y2(1,:).*tmp; 0107 0108 0109 % Rotate and translate back to original coordinates 0110 xo = tmp.*a+xo; 0111 yo = tmp.*b+yo; 0112 0113 % Mark special cases ................................... 0114 % First segments are degenerate to points 0115 if ~isempty(i_pt1) 0116 bool = ~s(i_pt1) & ~cc; 0117 i_pt12 = i_pt1(bool); 0118 i_pt1 = i_pt1(~bool); 0119 0120 bool = abs(tmp(i_pt1))<=tol; 0121 i_pt1_on = i_pt1(bool); 0122 i_pt1_off = i_pt1(~bool); 0123 0124 xo(i_pt1_on) = x1(1,i_pt1_on); 0125 yo(i_pt1_on) = y1(1,i_pt1_on); 0126 0127 oo = ones(size(i_pt1_off)); 0128 xo(i_pt1_off) = nan*oo; 0129 yo(i_pt1_off) = oo; 0130 end 0131 0132 % Second segments are degenerate to points ... 0133 if ~isempty(i_pt2) 0134 oo = ones(size(i_pt2_off)); 0135 xo(i_pt2_off) = nan*oo; 0136 yo(i_pt2_off) = 2*oo; 0137 end 0138 0139 % Both segments are degenerate ............... 0140 if ~isempty(i_pt12) 0141 bool = x1(i_pt12)==xo(i_pt12); 0142 i_pt12_on = i_pt12(bool); 0143 i_pt12_off = i_pt12(~bool); 0144 0145 xo(i_pt12_on) = x1(1,i_pt12_on); 0146 yo(i_pt12_on) = y1(1,i_pt12_on); 0147 0148 oo = ones(size(i_pt12_off)); 0149 xo(i_pt12_off) = nan*oo; 0150 yo(i_pt12_off) = 0*oo; 0151 end 0152 0153 % Parallel segments ......................... 0154 if ~isempty(i_par) 0155 oo = ones(size(i_par_on)); 0156 xo(i_par_on) = inf*oo; 0157 yo(i_par_on) = oo; 0158 0159 oo = ones(size(i_par_off)); 0160 xo(i_par_off) = inf*oo; 0161 yo(i_par_off) = 0*oo; 0162 end 0163 0164 0165 0166