polygonLoops

POLYGONLOOPS Divide a possibly self-intersecting polygon into a set of simple loops.

   LOOPS = polygonLoops(POLYGON);
   POLYGON is a polygone defined by a series of vertices,
   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:
       poly = [0 0;0 10;20 10;20 20;10 20;10 0];
       loops = polygonLoops(poly);
       figure(1); hold on;
       drawPolygon(loops);
       polygonArea(loops)

   See also
   polygons2d, polygonSelfIntersections

0001 function loops = polygonLoops(poly, varargin)
0002 %POLYGONLOOPS Divide a possibly self-intersecting polygon into a set of simple loops.
0003 %
0004 %   LOOPS = polygonLoops(POLYGON);
0005 %   POLYGON is a polygone defined by a series of vertices,
0006 %   LOOPS is a cell array of polygons, containing the same vertices of the
0007 %   original polygon, but no loop self-intersect, and no couple of loops
0008 %   intersect each other.
0009 %
0010 %   Example:
0011 %       poly = [0 0;0 10;20 10;20 20;10 20;10 0];
0012 %       loops = polygonLoops(poly);
0013 %       figure(1); hold on;
0014 %       drawPolygon(loops);
0015 %       polygonArea(loops)
0016 %
0017 %   See also
0018 %   polygons2d, polygonSelfIntersections
0019 %
0020 
0021 % ------
0022 % Author: David Legland
0023 % e-mail: david.legland@nantes.inra.fr
0024 % Created: 2009-06-15,    using Matlab 7.7.0.471 (R2008b)
0025 % Copyright 2009 INRA - Cepia Software Platform.
0026 
0027 
0028 % tolerance for detecting two vertices as equal
0029 tol = 1e-14;
0030 
0031 % parse optional arguments
0032 while length(varargin) > 1
0033     pname = varargin{1};
0034     if ~ischar(pname)
0035         error('Expect optional arguments as name-value pairs');
0036     end
0037     
0038     if strcmpi(pname, 'tolerance')
0039         tol = varargin{2};
0040     else
0041         error(['Unknown parameter name: ' pname]);
0042     end
0043     varargin(1:2) = [];
0044 end
0045 
0046 
0047 %% Initialisations
0048 
0049 % compute intersections
0050 [inters, pos1, pos2] = polygonSelfIntersections(poly, 'tolerance', tol);
0051 
0052 % case of a polygon without self-intersection
0053 if isempty(inters)
0054     loops = {poly};
0055     return;
0056 end
0057 
0058 % array for storing loops
0059 loops = cell(0, 1);
0060 
0061 % sort intersection points with respect to their position on the polygon
0062 [positions, order] = sortrows([pos1 pos2 ; pos2 pos1]);
0063 inters = [inters ; inters];
0064 inters = inters(order, :);
0065 
0066 
0067 %% First loop
0068 
0069 % initialize the beginning of the loop
0070 pos0 = 0;
0071 loop = polygonSubcurve(poly, pos0, positions(1, 1));
0072 loop(end, :) = inters(1,:);
0073 vertex = inters(1,:);
0074 
0075 % prepare iteration on positions
0076 pos = positions(1, 2);
0077 positions(1, :) = [];
0078 inters(1,:) = [];
0079 
0080 while true
0081     % index of next intersection point
0082     ind = find(positions(:,1) > pos, 1, 'first');
0083     
0084     % if not index is found, the current loop is complete
0085     if isempty(ind)
0086         break;
0087     end
0088 
0089     % compute the portion of curve between the two intersection points
0090     portion = polygonSubcurve(poly, pos, positions(ind, 1));
0091     
0092     % ensure extremities have been computed only once
0093     portion(1, :) = vertex;
0094     vertex = inters(ind, :);
0095     portion(end, :) = vertex;
0096     
0097     % add the current portion of curve
0098     loop = [loop; portion]; %#ok<AGROW>
0099     
0100     % update current position on the polygon
0101     pos = positions(ind, 2);
0102     
0103     % remove processed intersection
0104     positions(ind, :) = [];
0105     inters(ind,:) = [];
0106 end
0107 
0108 % append the last portion of curve
0109 loop = [loop ; polygonSubcurve(poly, pos, pos0)];
0110 
0111 % remove redundant vertices
0112 loop(sum(loop(1:end-1,:) == loop(2:end,:) ,2)==2, :) = [];
0113 if sum(diff(loop([1 end], :)) == 0) == 2
0114     loop(end, :) = [];
0115 end
0116 
0117 % add current loop to the list of loops
0118 loops{1} = loop;
0119 
0120 
0121 %% Other loops
0122 
0123 Nl = 1;
0124 while ~isempty(positions)
0125 
0126     % initialize the next loop
0127     loop    = [];
0128     pos0    = positions(1, 2);
0129     pos     = positions(1, 2);
0130     vertex  = inters(1,:);
0131     
0132     while true
0133         % index of next intersection point
0134         ind = find(positions(:,1) > pos, 1, 'first');
0135 
0136         % compute the portion of curve between the two intersection points
0137         portion = polygonSubcurve(poly, pos, positions(ind, 1));
0138         
0139         % ensure extremities have been computed only once
0140         portion(1, :) = vertex;
0141         vertex = inters(ind, :);
0142         portion(end, :) = vertex;
0143         
0144         % append the current portion of curve
0145         loop = [loop ; portion]; %#ok<AGROW>
0146 
0147         % update current position on the polygon
0148         pos = positions(ind, 2);
0149 
0150         % remove processed intersection
0151         positions(ind, :) = [];
0152         inters(ind,:) = [];
0153 
0154         % if not found, current loop is processed
0155         if pos == pos0
0156             break;
0157         end
0158     end
0159 
0160     % remove redundant vertices
0161     loop(sum(loop(1:end-1,:) == loop(2:end,:) ,2)==2, :) = []; %#ok<AGROW>
0162     if sum(diff(loop([1 end], :))==0) == 2
0163         loop(end, :) = []; 
0164     end
0165 
0166     % add current loop to the list of loops
0167     Nl = Nl + 1;
0168     loops{Nl} = loop;
0169 end