octave_packages/geometry-1.5.0/geom2d/intersectEdges.m

   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>
   4 %%
   5 %% All rights reserved.
   6 %% (simplified BSD License)
   7 %%
   8 %% Redistribution and use in source and binary forms, with or without
   9 %% modification, are permitted provided that the following conditions are met:
  10 %%
  11 %% 1. Redistributions of source code must retain the above copyright notice, this
  12 %%    list of conditions and the following disclaimer.
  13 %%
  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.
  17 %%
  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.
  29 %%
  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.
  33
  34 %% -*- texinfo -*-
  35 %% @deftypefn {Function File} {@var{point} = } intersectEdges (@var{edge1}, @var{edge2})
  36 %% Return all intersections between two set of edges
  37 %%
  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).
  42 %%
  43 %%   In case of colinear edges, returns [Inf Inf].
  44 %%   In case of parallel but not colinear edges, returns [NaN NaN].
  45 %%
  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
  49 %%   [N*2] array.
  50 %%
  51 %%   @seealso{edges2d, intersectLines}
  52 %% @end deftypefn
  53
  54 function point = intersectEdges(edge1, edge2)
  55   %% Initialisations
  56
  57   % ensure input arrays are same size
  58   N1  = size(edge1, 1);
  59   N2  = size(edge2, 1);
  60
  61   % ensure input have same size
  62   if N1~=N2
  63       if N1==1
  64           edge1 = repmat(edge1, [N2 1]);
  65           N1 = N2;
  66       elseif N2==1
  67           edge2 = repmat(edge2, [N1 1]);
  68       end
  69   end
  70
  71   % tolerance for precision
  72   tol = 1e-14;
  73
  74   % initialize result array
  75   x0  = zeros(N1, 1);
  76   y0  = zeros(N1, 1);
  77
  78
  79   %% Detect parallel and colinear cases
  80
  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));
  84
  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;
  89
  90   % Parallel edges have no intersection -> return [NaN NaN]
  91   x0(par & ~col) = NaN;
  92   y0(par & ~col) = NaN;
  93
  94
  95   %% Process colinear edges
  96
  97   % colinear edges may have 0, 1 or infinite intersection
  98   % Discrimnation based on position of edge2 vertices on edge1
  99   if sum(col)>0
 100       % array for storing results of colinear edges
 101       resCol = Inf*ones(size(col));
 102
 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, :));
 106
 107       % control location of vertices: we want t1<t2
 108       if t1>t2
 109           tmp = t1;
 110           t1  = t2;
 111           t2  = tmp;
 112       end
 113
 114       % edge totally before first vertex or totally after last vertex
 115       resCol(col(t2<-eps))  = NaN;
 116       resCol(col(t1>1+eps)) = NaN;
 117
 118       % set up result into point coordinate
 119       x0(col) = resCol;
 120       y0(col) = resCol;
 121
 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);
 126
 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);
 131   end
 132
 133
 134   %% Process non parallel cases
 135
 136   % process intersecting edges whose interecting lines intersect
 137   i = find(~par);
 138
 139   % use a test to avoid process empty arrays
 140   if sum(i)>0
 141       % extract base parameters of supporting lines for non-parallel edges
 142       x1  = edge1(i,1);
 143       y1  = edge1(i,2);
 144       dx1 = edge1(i,3)-x1;
 145       dy1 = edge1(i,4)-y1;
 146       x2  = edge2(i,1);
 147       y2  = edge2(i,2);
 148       dx2 = edge2(i,3)-x2;
 149       dy2 = edge2(i,4)-y2;
 150
 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;
 155
 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);
 160
 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;
 165       x0(i(out)) = NaN;
 166       y0(i(out)) = NaN;
 167   end
 168
 169
 170   %% format output arguments
 171
 172   point = [x0 y0];
 173
 174 endfunction
 175