1 %% Copyright (c) 2011, INRA
2 %% 2003-2011, David Legland <david.legland@grignon.inra.fr>
3 %% 2011 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
5 %% All rights reserved.
6 %% (simplified BSD License)
8 %% Redistribution and use in source and binary forms, with or without
9 %% modification, are permitted provided that the following conditions are met:
11 %% 1. Redistributions of source code must retain the above copyright notice, this
12 %% list of conditions and the following disclaimer.
14 %% 2. Redistributions in binary form must reproduce the above copyright notice,
15 %% this list of conditions and the following disclaimer in the documentation
16 %% and/or other materials provided with the distribution.
18 %% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 %% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 %% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21 %% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
22 %% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
23 %% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
24 %% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25 %% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
26 %% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
27 %% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
28 %% POSSIBILITY OF SUCH DAMAGE.
30 %% The views and conclusions contained in the software and documentation are
31 %% those of the authors and should not be interpreted as representing official
32 %% policies, either expressed or implied, of copyright holder.
35 %% @deftypefn {Function File} {@var{point} = } intersectEdges (@var{edge1}, @var{edge2})
36 %% Return all intersections between two set of edges
38 %% P = intersectEdges(E1, E2);
39 %% returns the intersection point of lines L1 and L2. E1 and E2 are 1-by-4
40 %% arrays, containing parametric representation of each edge (in the form
41 %% [x1 y1 x2 y2], see 'createEdge' for details).
43 %% In case of colinear edges, returns [Inf Inf].
44 %% In case of parallel but not colinear edges, returns [NaN NaN].
46 %% If each input is [N*4] array, the result is a [N*2] array containing
47 %% intersections of each couple of edges.
48 %% If one of the input has N rows and the other 1 row, the result is a
51 %% @seealso{edges2d, intersectLines}
54 function point = intersectEdges(edge1, edge2)
57 % ensure input arrays are same size
61 % ensure input have same size
64 edge1 = repmat(edge1, [N2 1]);
67 edge2 = repmat(edge2, [N1 1]);
71 % tolerance for precision
74 % initialize result array
79 %% Detect parallel and colinear cases
81 % indices of parallel edges
82 %par = abs(dx1.*dy2 - dx1.*dy2)<tol;
83 par = isParallel(edge1(:,3:4)-edge1(:,1:2), edge2(:,3:4)-edge2(:,1:2));
85 % indices of colinear edges
86 % equivalent to: |(x2-x1)*dy1 - (y2-y1)*dx1| < eps
87 col = abs( (edge2(:,1)-edge1(:,1)).*(edge1(:,4)-edge1(:,2)) - ...
88 (edge2(:,2)-edge1(:,2)).*(edge1(:,3)-edge1(:,1)) )<tol & par;
90 % Parallel edges have no intersection -> return [NaN NaN]
95 %% Process colinear edges
97 % colinear edges may have 0, 1 or infinite intersection
98 % Discrimnation based on position of edge2 vertices on edge1
100 % array for storing results of colinear edges
101 resCol = Inf*ones(size(col));
103 % compute position of edge2 vertices wrt edge1
104 t1 = edgePosition(edge2(col, 1:2), edge1(col, :));
105 t2 = edgePosition(edge2(col, 3:4), edge1(col, :));
107 % control location of vertices: we want t1<t2
114 % edge totally before first vertex or totally after last vertex
115 resCol(col(t2<-eps)) = NaN;
116 resCol(col(t1>1+eps)) = NaN;
118 % set up result into point coordinate
122 % touches on first point of first edge
123 touch = col(abs(t2)<1e-14);
124 x0(touch) = edge1(touch, 1);
125 y0(touch) = edge1(touch, 2);
127 % touches on second point of first edge
128 touch = col(abs(t1-1)<1e-14);
129 x0(touch) = edge1(touch, 3);
130 y0(touch) = edge1(touch, 4);
134 %% Process non parallel cases
136 % process intersecting edges whose interecting lines intersect
139 % use a test to avoid process empty arrays
141 % extract base parameters of supporting lines for non-parallel edges
151 % compute intersection points of supporting lines
152 delta = (dx2.*dy1-dx1.*dy2);
153 x0(i) = ((y2-y1).*dx1.*dx2 + x1.*dy1.*dx2 - x2.*dy2.*dx1) ./ delta;
154 y0(i) = ((x2-x1).*dy1.*dy2 + y1.*dx1.*dy2 - y2.*dx2.*dy1) ./ -delta;
156 % compute position of intersection points on each edge
157 % t1 is position on edge1, t2 is position on edge2
158 t1 = ((y0(i)-y1).*dy1 + (x0(i)-x1).*dx1) ./ (dx1.*dx1+dy1.*dy1);
159 t2 = ((y0(i)-y2).*dy2 + (x0(i)-x2).*dx2) ./ (dx2.*dx2+dy2.*dy2);
161 % check position of points on edges.
162 % it should be comprised between 0 and 1 for both t1 and t2.
163 % if not, the edges do not intersect
164 out = t1<-tol | t1>1+tol | t2<-tol | t2>1+tol;
170 %% format output arguments