INTERSECTEDGES Return all intersections between two set of edges. P = intersectEdges(E1, E2); returns the intersection point of edges E1 and E2. E1 and E2 are 1-by-4 arrays, containing parametric representation of each edge (in the form [x1 y1 x2 y2], see 'createEdge' for details). In case of colinear edges, the result P contains [Inf Inf]. In case of parallel but not colinear edges, the result P contains [NaN NaN]. If each input is N-by-4 array, the result is a N-by-2 array containing the intersection of each couple of edges. If one of the input has N rows and the other 1 row, the result is a N-by-2 array. P = intersectEdges(E1, E2, TOL); Specifies a tolerance parameter to determine parallel and colinear edges, and if a point belongs to an edge or not. The latter test is performed on the relative position of the intersection point over the edge, that should lie within [-TOL; 1+TOL]. See also edges2d, intersectLines
0001 function point = intersectEdges(edge1, edge2, varargin) 0002 %INTERSECTEDGES Return all intersections between two set of edges. 0003 % 0004 % P = intersectEdges(E1, E2); 0005 % returns the intersection point of edges E1 and E2. 0006 % E1 and E2 are 1-by-4 arrays, containing parametric representation of 0007 % each edge (in the form [x1 y1 x2 y2], see 'createEdge' for details). 0008 % 0009 % In case of colinear edges, the result P contains [Inf Inf]. 0010 % In case of parallel but not colinear edges, the result P contains 0011 % [NaN NaN]. 0012 % 0013 % If each input is N-by-4 array, the result is a N-by-2 array containing 0014 % the intersection of each couple of edges. 0015 % If one of the input has N rows and the other 1 row, the result is a 0016 % N-by-2 array. 0017 % 0018 % P = intersectEdges(E1, E2, TOL); 0019 % Specifies a tolerance parameter to determine parallel and colinear 0020 % edges, and if a point belongs to an edge or not. The latter test is 0021 % performed on the relative position of the intersection point over the 0022 % edge, that should lie within [-TOL; 1+TOL]. 0023 % 0024 % See also 0025 % edges2d, intersectLines 0026 % 0027 0028 % ------ 0029 % Author: David Legland 0030 % E-mail: david.legland@inrae.fr 0031 % Created: 2003-10-31 0032 % Copyright 2003-2024 INRA - Cepia Software Platform 0033 0034 % tolerance for precision 0035 tol = 1e-14; 0036 0037 if nargin > 2 0038 tol = varargin{1}; 0039 end 0040 0041 0042 %% Initialisations 0043 0044 % ensure input arrays are same size 0045 N1 = size(edge1, 1); 0046 N2 = size(edge2, 1); 0047 0048 % ensure input have same size 0049 if N1 ~= N2 0050 if N1 == 1 0051 edge1 = repmat(edge1, [N2 1]); 0052 N1 = N2; 0053 elseif N2 == 1 0054 edge2 = repmat(edge2, [N1 1]); 0055 end 0056 end 0057 0058 % initialize result array 0059 x0 = zeros(N1, 1); 0060 y0 = zeros(N1, 1); 0061 0062 0063 %% Detect parallel and colinear cases 0064 0065 % indices of parallel edges 0066 %par = abs(dx1.*dy2 - dx1.*dy2)<tol; 0067 par = isParallel(edge1(:,3:4)-edge1(:,1:2), edge2(:,3:4)-edge2(:,1:2)); 0068 0069 % indices of colinear edges 0070 % equivalent to: |(x2-x1)*dy1 - (y2-y1)*dx1| < eps 0071 col = abs( (edge2(:,1)-edge1(:,1)) .* (edge1(:,4)-edge1(:,2)) - ... 0072 (edge2(:,2)-edge1(:,2)) .* (edge1(:,3)-edge1(:,1)) ) < tol & par; 0073 0074 % Parallel edges have no intersection -> return [NaN NaN] 0075 x0(par & ~col) = NaN; 0076 y0(par & ~col) = NaN; 0077 0078 0079 %% Process colinear edges 0080 0081 % colinear edges may have 0, 1 or infinite intersection 0082 % Setup result intersection point as follow: 0083 % * no intersection -> [NaN NaN] 0084 % * partial overlap -> [Inf Inf] 0085 % * touches at extremity -> extremity coordinates 0086 % Discrimnation based on position of edge2 vertices on edge1 0087 if sum(col) > 0 0088 % array for storing results of colinear edges 0089 resCol = Inf * ones(size(col)); 0090 colInds = find(col); 0091 0092 % compute position of edge2 vertices wrt edge1 0093 t1 = edgePosition(edge2(col, 1:2), edge1(col, :), 'diag'); 0094 t2 = edgePosition(edge2(col, 3:4), edge1(col, :), 'diag'); 0095 0096 % control location of vertices: we want t1<t2 0097 swap = t1 > t2; 0098 tmp = t1(swap); 0099 t1(swap) = t2(swap); 0100 t2(swap) = tmp; 0101 0102 % edge totally before first vertex or totally after last vertex 0103 resCol(colInds(t2 < -tol)) = NaN; 0104 resCol(colInds(t1 > 1+tol)) = NaN; 0105 0106 % set up result into point coordinate 0107 x0(col) = resCol(col); 0108 y0(col) = resCol(col); 0109 0110 % touches on first point of first edge 0111 touch = colInds(abs(t2) < tol); 0112 x0(touch) = edge1(touch, 1); 0113 y0(touch) = edge1(touch, 2); 0114 0115 % touches on second point of first edge 0116 touch = colInds(abs(t1-1) < tol); 0117 x0(touch) = edge1(touch, 3); 0118 y0(touch) = edge1(touch, 4); 0119 end 0120 0121 0122 %% Process non parallel cases 0123 0124 % process edges whose supporting lines intersect 0125 i = find(~par); 0126 0127 % use a test to avoid process empty arrays 0128 if sum(i) > 0 0129 % extract base parameters of supporting lines for non-parallel edges 0130 x1 = edge1(i,1); 0131 y1 = edge1(i,2); 0132 dx1 = edge1(i,3) - x1; 0133 dy1 = edge1(i,4) - y1; 0134 x2 = edge2(i,1); 0135 y2 = edge2(i,2); 0136 dx2 = edge2(i,3) - x2; 0137 dy2 = edge2(i,4) - y2; 0138 0139 % compute intersection points of supporting lines 0140 delta = (dx2.*dy1 - dx1.*dy2); 0141 x0(i) = ((y2-y1).*dx1.*dx2 + x1.*dy1.*dx2 - x2.*dy2.*dx1) ./ delta; 0142 y0(i) = ((x2-x1).*dy1.*dy2 + y1.*dx1.*dy2 - y2.*dx2.*dy1) ./ -delta; 0143 0144 % compute position of intersection points on each edge 0145 % t1 is position on edge1, t2 is position on edge2 0146 t1 = ((y0(i)-y1).*dy1 + (x0(i)-x1).*dx1) ./ (dx1.*dx1+dy1.*dy1); 0147 t2 = ((y0(i)-y2).*dy2 + (x0(i)-x2).*dx2) ./ (dx2.*dx2+dy2.*dy2); 0148 0149 % check position of points on edges. 0150 % it should be comprised between 0 and 1 for both t1 and t2. 0151 % if not, the edges do not intersect 0152 out = t1<-tol | t1>1+tol | t2<-tol | t2>1+tol; 0153 x0(i(out)) = NaN; 0154 y0(i(out)) = NaN; 0155 end 0156 0157 0158 %% format output arguments 0159 0160 point = [x0 y0]; 0161