Add a useful package (from Source forge) for octave

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

diff --git a/octave_packages/geometry-1.5.0/polygons2d/polygonLoops.m b/octave_packages/geometry-1.5.0/polygons2d/polygonLoops.m
new file mode 100644 (file)
index 0000000..e024f6b
--- /dev/null
+++ b/octave_packages/geometry-1.5.0/polygons2d/polygonLoops.m
@@ -0,0 +1,143 @@
+## 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{loops} = } polygonLoops (@var{poly})
+## Divide a possibly self-intersecting polygon into a set of simple loops
+##
+##   @var{poly} is a polygone defined by a series of vertices,
+##   @var{loops} is a cell array of polygons, containing the same vertices of the
+##   original polygon, but no loop self-intersect, and no couple of loops
+##   intersect each other.
+##
+##   Example:
+## @example
+##       poly = [0 0;0 10;20 10;20 20;10 20;10 0];
+##       loops = polygonLoops(poly);
+##       figure(1); hold on;
+##       drawPolygon(loops);
+##       polygonArea(loops)
+## @end example
+##
+## @seealso{polygons2d, polygonSelfIntersections}
+## @end deftypefn
+
+function loops = polygonLoops(poly)
+
+  ## Initialisations
+
+  # compute intersections
+  [inters pos1 pos2] = polygonSelfIntersections(poly);
+
+  # case of a polygon without self-intersection
+  if isempty(inters)
+      loops = {poly};
+      return;
+  end
+
+  # array for storing loops
+  loops = cell(0, 1);
+
+  positions = sortrows([pos1 pos2;pos2 pos1]);
+
+
+  ## First loop
+
+  pos0 = 0;
+  loop = polygonSubcurve(poly, pos0, positions(1, 1));
+  pos = positions(1, 2);
+  positions(1, :) = [];
+
+  while true
+      # index of next intersection point
+      ind = find(positions(:,1)>pos, 1, 'first');
+      
+      # if not found, break
+      if isempty(ind)
+          break;
+      end
+      
+      # add portion of curve
+      loop = [loop;polygonSubcurve(poly, pos, positions(ind, 1))]; ##ok<AGROW>
+      
+      # look for next intersection point
+      pos = positions(ind, 2);
+      positions(ind, :) = [];
+  end
+
+  # add the last portion of curve
+  loop = [loop;polygonSubcurve(poly, pos, pos0)];
+
+  # remove redundant vertices
+  loop(sum(loop(1:end-1,:) == loop(2:end,:) ,2)==2, :) = [];
+  if sum(diff(loop([1 end], :))==0)==2
+      loop(end, :) = [];
+  end
+
+  # add current loop to the list of loops
+  loops{1} = loop;
+
+
+  ## Other loops
+
+  Nl = 1;
+  while ~isempty(positions)
+
+      loop    = [];
+      pos0    = positions(1, 2);
+      pos     = positions(1, 2);
+      
+      while true
+          # index of next intersection point
+          ind = find(positions(:,1)>pos, 1, 'first');
+
+          # add portion of curve
+          loop = [loop;polygonSubcurve(poly, pos, positions(ind, 1))]; ##ok<AGROW>
+
+          # look for next intersection point
+          pos = positions(ind, 2);
+          positions(ind, :) = [];
+
+          # if not found, break
+          if pos==pos0
+              break;
+          end
+      end
+
+      # remove redundant vertices
+      loop(sum(loop(1:end-1,:) == loop(2:end,:) ,2)==2, :) = []; ##ok<AGROW>
+      if sum(diff(loop([1 end], :))==0)==2
+          loop(end, :) = []; ##ok<AGROW>
+      end
+
+      # add current loop to the list of loops
+      Nl = Nl + 1;
+      loops{Nl} = loop;
+  end
+
+endfunction