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@nantes.inra.fr 0031 % INRA - Cepia Software Platform 0032 % created the 31/10/2003. 0033 % 0034 0035 % HISTORY 0036 % 19/02/2004 add support for multiple edges. 0037 % 15/08/2005 rewrite a lot, and create unit tests 0038 % 08/03/2007 update doc 0039 % 21/09/2009 fix bug for edge arrays (thanks to Miquel Cubells) 0040 % 03/08/2010 fix another bug for edge arrays (thanks to Reto Zingg) 0041 % 20/02/2013 fix bug for management of parallel edges (thanks to Nicolaj 0042 % Kirchhof) 0043 % 19/07/2016 add support for tolerance 0044 0045 % tolerance for precision 0046 tol = 1e-14; 0047 0048 if nargin > 2 0049 tol = varargin{1}; 0050 end 0051 0052 0053 %% Initialisations 0054 0055 % ensure input arrays are same size 0056 N1 = size(edge1, 1); 0057 N2 = size(edge2, 1); 0058 0059 % ensure input have same size 0060 if N1 ~= N2 0061 if N1 == 1 0062 edge1 = repmat(edge1, [N2 1]); 0063 N1 = N2; 0064 elseif N2 == 1 0065 edge2 = repmat(edge2, [N1 1]); 0066 end 0067 end 0068 0069 % initialize result array 0070 x0 = zeros(N1, 1); 0071 y0 = zeros(N1, 1); 0072 0073 0074 %% Detect parallel and colinear cases 0075 0076 % indices of parallel edges 0077 %par = abs(dx1.*dy2 - dx1.*dy2)<tol; 0078 par = isParallel(edge1(:,3:4)-edge1(:,1:2), edge2(:,3:4)-edge2(:,1:2)); 0079 0080 % indices of colinear edges 0081 % equivalent to: |(x2-x1)*dy1 - (y2-y1)*dx1| < eps 0082 col = abs( (edge2(:,1)-edge1(:,1)) .* (edge1(:,4)-edge1(:,2)) - ... 0083 (edge2(:,2)-edge1(:,2)) .* (edge1(:,3)-edge1(:,1)) ) < tol & par; 0084 0085 % Parallel edges have no intersection -> return [NaN NaN] 0086 x0(par & ~col) = NaN; 0087 y0(par & ~col) = NaN; 0088 0089 0090 %% Process colinear edges 0091 0092 % colinear edges may have 0, 1 or infinite intersection 0093 % Discrimnation based on position of edge2 vertices on edge1 0094 if sum(col) > 0 0095 % array for storing results of colinear edges 0096 resCol = Inf * ones(size(col)); 0097 0098 % compute position of edge2 vertices wrt edge1 0099 t1 = edgePosition(edge2(col, 1:2), edge1(col, :)); 0100 t2 = edgePosition(edge2(col, 3:4), edge1(col, :)); 0101 0102 % control location of vertices: we want t1<t2 0103 if t1 > t2 0104 tmp = t1; 0105 t1 = t2; 0106 t2 = tmp; 0107 end 0108 0109 % edge totally before first vertex or totally after last vertex 0110 resCol(col(t2 < -tol)) = NaN; 0111 resCol(col(t1 > 1+tol)) = NaN; 0112 0113 % set up result into point coordinate 0114 x0(col) = resCol(col); 0115 y0(col) = resCol(col); 0116 0117 % touches on first point of first edge 0118 touch = col(abs(t2) < tol); 0119 x0(touch) = edge1(touch, 1); 0120 y0(touch) = edge1(touch, 2); 0121 0122 % touches on second point of first edge 0123 touch = col(abs(t1-1) < tol); 0124 x0(touch) = edge1(touch, 3); 0125 y0(touch) = edge1(touch, 4); 0126 end 0127 0128 0129 %% Process non parallel cases 0130 0131 % process edges whose supporting lines intersect 0132 i = find(~par); 0133 0134 % use a test to avoid process empty arrays 0135 if sum(i) > 0 0136 % extract base parameters of supporting lines for non-parallel edges 0137 x1 = edge1(i,1); 0138 y1 = edge1(i,2); 0139 dx1 = edge1(i,3) - x1; 0140 dy1 = edge1(i,4) - y1; 0141 x2 = edge2(i,1); 0142 y2 = edge2(i,2); 0143 dx2 = edge2(i,3) - x2; 0144 dy2 = edge2(i,4) - y2; 0145 0146 % compute intersection points of supporting lines 0147 delta = (dx2.*dy1 - dx1.*dy2); 0148 x0(i) = ((y2-y1).*dx1.*dx2 + x1.*dy1.*dx2 - x2.*dy2.*dx1) ./ delta; 0149 y0(i) = ((x2-x1).*dy1.*dy2 + y1.*dx1.*dy2 - y2.*dx2.*dy1) ./ -delta; 0150 0151 % compute position of intersection points on each edge 0152 % t1 is position on edge1, t2 is position on edge2 0153 t1 = ((y0(i)-y1).*dy1 + (x0(i)-x1).*dx1) ./ (dx1.*dx1+dy1.*dy1); 0154 t2 = ((y0(i)-y2).*dy2 + (x0(i)-x2).*dx2) ./ (dx2.*dx2+dy2.*dy2); 0155 0156 % check position of points on edges. 0157 % it should be comprised between 0 and 1 for both t1 and t2. 0158 % if not, the edges do not intersect 0159 out = t1<-tol | t1>1+tol | t2<-tol | t2>1+tol; 0160 x0(i(out)) = NaN; 0161 y0(i(out)) = NaN; 0162 end 0163 0164 0165 %% format output arguments 0166 0167 point = [x0 y0]; 0168