octave_packages/geometry-1.5.0/polygons2d/polygonLoops.m

   1 ## Copyright (C) 2003-2011 David Legland <david.legland@grignon.inra.fr>
   2 ## Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
   3 ## All rights reserved.
   4 ##
   5 ## Redistribution and use in source and binary forms, with or without
   6 ## modification, are permitted provided that the following conditions are met:
   7 ##
   8 ##     1 Redistributions of source code must retain the above copyright notice,
   9 ##       this list of conditions and the following disclaimer.
  10 ##     2 Redistributions in binary form must reproduce the above copyright
  11 ##       notice, this list of conditions and the following disclaimer in the
  12 ##       documentation and/or other materials provided with the distribution.
  13 ##
  14 ## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
  15 ## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  16 ## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  17 ## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
  18 ## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  19 ## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
  20 ## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  21 ## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
  22 ## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  23 ## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  24 ##
  25 ## The views and conclusions contained in the software and documentation are
  26 ## those of the authors and should not be interpreted as representing official
  27 ## policies, either expressed or implied, of the copyright holders.
  28
  29 ## -*- texinfo -*-
  30 ## @deftypefn {Function File} {@var{loops} = } polygonLoops (@var{poly})
  31 ## Divide a possibly self-intersecting polygon into a set of simple loops
  32 ##
  33 ##   @var{poly} is a polygone defined by a series of vertices,
  34 ##   @var{loops} is a cell array of polygons, containing the same vertices of the
  35 ##   original polygon, but no loop self-intersect, and no couple of loops
  36 ##   intersect each other.
  37 ##
  38 ##   Example:
  39 ## @example
  40 ##       poly = [0 0;0 10;20 10;20 20;10 20;10 0];
  41 ##       loops = polygonLoops(poly);
  42 ##       figure(1); hold on;
  43 ##       drawPolygon(loops);
  44 ##       polygonArea(loops)
  45 ## @end example
  46 ##
  47 ## @seealso{polygons2d, polygonSelfIntersections}
  48 ## @end deftypefn
  49
  50 function loops = polygonLoops(poly)
  51
  52   ## Initialisations
  53
  54   # compute intersections
  55   [inters pos1 pos2] = polygonSelfIntersections(poly);
  56
  57   # case of a polygon without self-intersection
  58   if isempty(inters)
  59       loops = {poly};
  60       return;
  61   end
  62
  63   # array for storing loops
  64   loops = cell(0, 1);
  65
  66   positions = sortrows([pos1 pos2;pos2 pos1]);
  67
  68
  69   ## First loop
  70
  71   pos0 = 0;
  72   loop = polygonSubcurve(poly, pos0, positions(1, 1));
  73   pos = positions(1, 2);
  74   positions(1, :) = [];
  75
  76   while true
  77       # index of next intersection point
  78       ind = find(positions(:,1)>pos, 1, 'first');
  79
  80       # if not found, break
  81       if isempty(ind)
  82           break;
  83       end
  84
  85       # add portion of curve
  86       loop = [loop;polygonSubcurve(poly, pos, positions(ind, 1))]; ##ok<AGROW>
  87
  88       # look for next intersection point
  89       pos = positions(ind, 2);
  90       positions(ind, :) = [];
  91   end
  92
  93   # add the last portion of curve
  94   loop = [loop;polygonSubcurve(poly, pos, pos0)];
  95
  96   # remove redundant vertices
  97   loop(sum(loop(1:end-1,:) == loop(2:end,:) ,2)==2, :) = [];
  98   if sum(diff(loop([1 end], :))==0)==2
  99       loop(end, :) = [];
 100   end
 101
 102   # add current loop to the list of loops
 103   loops{1} = loop;
 104
 105
 106   ## Other loops
 107
 108   Nl = 1;
 109   while ~isempty(positions)
 110
 111       loop    = [];
 112       pos0    = positions(1, 2);
 113       pos     = positions(1, 2);
 114
 115       while true
 116           # index of next intersection point
 117           ind = find(positions(:,1)>pos, 1, 'first');
 118
 119           # add portion of curve
 120           loop = [loop;polygonSubcurve(poly, pos, positions(ind, 1))]; ##ok<AGROW>
 121
 122           # look for next intersection point
 123           pos = positions(ind, 2);
 124           positions(ind, :) = [];
 125
 126           # if not found, break
 127           if pos==pos0
 128               break;
 129           end
 130       end
 131
 132       # remove redundant vertices
 133       loop(sum(loop(1:end-1,:) == loop(2:end,:) ,2)==2, :) = []; ##ok<AGROW>
 134       if sum(diff(loop([1 end], :))==0)==2
 135           loop(end, :) = []; ##ok<AGROW>
 136       end
 137
 138       # add current loop to the list of loops
 139       Nl = Nl + 1;
 140       loops{Nl} = loop;
 141   end
 142
 143 endfunction