Add a useful package (from Source forge) for octave

[CreaPhase.git] / octave_packages / geometry-1.5.0 / polygons2d / polylineSelfIntersections.m

diff --git a/octave_packages/geometry-1.5.0/polygons2d/polylineSelfIntersections.m b/octave_packages/geometry-1.5.0/polygons2d/polylineSelfIntersections.m
new file mode 100644 (file)
index 0000000..d759be9
--- /dev/null
+++ b/octave_packages/geometry-1.5.0/polygons2d/polylineSelfIntersections.m
@@ -0,0 +1,153 @@
+## Copyright (C) 2003-2011 David Legland <david.legland@grignon.inra.fr>
+## Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
+## All rights reserved.
+##
+## Redistribution and use in source and binary forms, with or without
+## modification, are permitted provided that the following conditions are met:
+##
+##     1 Redistributions of source code must retain the above copyright notice,
+##       this list of conditions and the following disclaimer.
+##     2 Redistributions in binary form must reproduce the above copyright
+##       notice, this list of conditions and the following disclaimer in the
+##       documentation and/or other materials provided with the distribution.
+##
+## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
+## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
+## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+##
+## The views and conclusions contained in the software and documentation are
+## those of the authors and should not be interpreted as representing official
+## policies, either expressed or implied, of the copyright holders.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{pts} = } polylineSelfIntersections (@var{poly})
+##   Find self-intersections points of a polyline
+##
+##   Return the position of self intersection points
+##
+##   Also return the 2 positions of each intersection point (the position
+##   when meeting point for first time, then position when meeting point
+##   for the second time).
+##
+##   Example
+## @example
+##       # use a gamma-shaped polyline
+##       poly = [0 0;0 10;20 10;20 20;10 20;10 0];
+##       polylineSelfIntersections(poly)
+##       ans = 
+##           10 10
+##
+##       # use a 'S'-shaped polyline
+##       poly = [10 0;0 0;0 10;20 10;20 20;10 20];
+##       polylineSelfIntersections(poly)
+##       ans = 
+##           10 10
+## @end example
+##   
+##  @seealso{polygons2d, intersectPolylines, polygonSelfIntersections}
+## @end deftypefn
+
+function varargout = polylineSelfIntersections(poly, varargin)
+
+  ## Initialisations
+
+  # determine whether the polyline is closed
+  closed = false;
+  if ~isempty(varargin)
+      closed = varargin{1};
+      if ischar(closed)
+          if strcmp(closed, 'closed')
+              closed = true;
+          elseif strcmp(closed, 'open')
+              closed = false;
+          end
+      end
+  end
+
+  # if polyline is closed, ensure the last point equals the first one
+  if closed
+      if sum(abs(poly(end, :)-poly(1,:))<1e-14)~=2
+          poly = [poly; poly(1,:)];
+      end
+  end
+
+  # arrays for storing results
+  points  = zeros(0, 2);
+  pos1    = zeros(0, 1);
+  pos2    = zeros(0, 1);
+
+  # number of vertices
+  Nv = size(poly, 1);
+
+
+  ## Main processing
+
+  # index of current intersection
+  ip = 0;
+
+  # iterate over each couple of edge ( (N-1)*(N-2)/2 iterations)
+  for i=1:Nv-2
+      # create first edge
+      edge1 = [poly(i, :) poly(i+1, :)];
+      for j=i+2:Nv-1
+          # create second edge
+          edge2 = [poly(j, :) poly(j+1, :)];
+
+          # check conditions on bounding boxes
+          if min(edge1([1 3]))>max(edge2([1 3]))
+              continue;
+          end
+          if max(edge1([1 3]))<min(edge2([1 3]))
+              continue;
+          end
+          if min(edge1([2 4]))>max(edge2([2 4]))
+              continue;
+          end
+          if max(edge1([2 4]))<min(edge2([2 4]))
+              continue;
+          end
+          
+          # compute intersection point
+          inter = intersectEdges(edge1, edge2);
+          
+          if sum(isfinite(inter))==2
+              # add point to the list
+              ip = ip + 1;
+              points(ip, :) = inter;
+              
+              # also compute positions on the polyline
+              pos1(ip, 1) = i+edgePosition(inter, edge1)-1;
+              pos2(ip, 1) = j+edgePosition(inter, edge2)-1;
+          end
+      end
+  end
+
+  # if polyline is closed, the first vertex was found as an intersection, so
+  # we need to remove it
+  if closed
+      dist = distancePoints(points, poly(1,:));
+      [minDist ind] = min(dist); ##ok<ASGLU>
+      points(ind,:) = [];
+      pos1(ind)   = [];
+      pos2(ind)   = [];
+  end
+
+  ## Post-processing
+
+  # process output arguments
+  if nargout<=1
+      varargout{1} = points;
+  elseif nargout==3
+      varargout{1} = points;
+      varargout{2} = pos1;
+      varargout{3} = pos2;
+  end
+
+endfunction