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@inrae.fr
0024 % Created: 2009-06-15, using Matlab 7.7.0.471 (R2008b)
0025 % Copyright 2009-2024 INRA - Cepia Software Platform
0026 
0027 % tolerance for detecting two vertices as equal
0028 tol = 1e-14;
0029 
0030 % parse optional arguments
0031 while length(varargin) > 1
0032     pname = varargin{1};
0033     if ~ischar(pname)
0034         error('Expect optional arguments as name-value pairs');
0035     end
0036     
0037     if strcmpi(pname, 'tolerance')
0038         tol = varargin{2};
0039     else
0040         error(['Unknown parameter name: ' pname]);
0041     end
0042     varargin(1:2) = [];
0043 end
0044 
0045 
0046 %% Initialisations
0047 
0048 % compute intersections
0049 [inters, pos1, pos2] = polygonSelfIntersections(poly, 'tolerance', tol);
0050 
0051 % case of a polygon without self-intersection
0052 if isempty(inters)
0053     loops = {poly};
0054     return;
0055 end
0056 
0057 % array for storing loops
0058 loops = cell(0, 1);
0059 
0060 % sort intersection points with respect to their position on the polygon
0061 [positions, order] = sortrows([pos1 pos2 ; pos2 pos1]);
0062 inters = [inters ; inters];
0063 inters = inters(order, :);
0064 
0065 
0066 %% First loop
0067 
0068 % initialize the beginning of the loop
0069 pos0 = 0;
0070 loop = polygonSubcurve(poly, pos0, positions(1, 1));
0071 loop(end, :) = inters(1,:);
0072 vertex = inters(1,:);
0073 
0074 % prepare iteration on positions
0075 pos = positions(1, 2);
0076 positions(1, :) = [];
0077 inters(1,:) = [];
0078 
0079 while true
0080     % index of next intersection point
0081     ind = find(positions(:,1) > pos, 1, 'first');
0082     
0083     % if not index is found, the current loop is complete
0084     if isempty(ind)
0085         break;
0086     end
0087 
0088     % compute the portion of curve between the two intersection points
0089     portion = polygonSubcurve(poly, pos, positions(ind, 1));
0090     
0091     % ensure extremities have been computed only once
0092     portion(1, :) = vertex;
0093     vertex = inters(ind, :);
0094     portion(end, :) = vertex;
0095     
0096     % add the current portion of curve
0097     loop = [loop; portion]; %#ok<AGROW>
0098     
0099     % update current position on the polygon
0100     pos = positions(ind, 2);
0101     
0102     % remove processed intersection
0103     positions(ind, :) = [];
0104     inters(ind,:) = [];
0105 end
0106 
0107 % append the last portion of curve
0108 loop = [loop ; polygonSubcurve(poly, pos, pos0)];
0109 
0110 % remove redundant vertices
0111 loop(sum(loop(1:end-1,:) == loop(2:end,:) ,2)==2, :) = [];
0112 if sum(diff(loop([1 end], :)) == 0) == 2
0113     loop(end, :) = [];
0114 end
0115 
0116 % add current loop to the list of loops
0117 loops{1} = loop;
0118 
0119 
0120 %% Other loops
0121 
0122 Nl = 1;
0123 while ~isempty(positions)
0124 
0125     % initialize the next loop
0126     loop    = [];
0127     pos0    = positions(1, 2);
0128     pos     = positions(1, 2);
0129     vertex  = inters(1,:);
0130     
0131     while true
0132         % index of next intersection point
0133         ind = find(positions(:,1) > pos, 1, 'first');
0134 
0135         % compute the portion of curve between the two intersection points
0136         portion = polygonSubcurve(poly, pos, positions(ind, 1));
0137         
0138         % ensure extremities have been computed only once
0139         portion(1, :) = vertex;
0140         vertex = inters(ind, :);
0141         portion(end, :) = vertex;
0142         
0143         % append the current portion of curve
0144         loop = [loop ; portion]; %#ok<AGROW>
0145 
0146         % update current position on the polygon
0147         pos = positions(ind, 2);
0148 
0149         % remove processed intersection
0150         positions(ind, :) = [];
0151         inters(ind,:) = [];
0152 
0153         % if not found, current loop is processed
0154         if pos == pos0
0155             break;
0156         end
0157     end
0158 
0159     % remove redundant vertices
0160     loop(sum(loop(1:end-1,:) == loop(2:end,:) ,2)==2, :) = []; %#ok<AGROW>
0161     if sum(diff(loop([1 end], :))==0) == 2
0162         loop(end, :) = []; 
0163     end
0164 
0165     % add current loop to the list of loops
0166     Nl = Nl + 1;
0167     loops{Nl} = loop;
0168 end