POLYLINESELFINTERSECTIONS Find self-intersection points of a polyline.
Computes self-intersections of a polyline, eventually specifying if
polyline is closed or open, and eventually returning position of
intersection points on polyline.
For common use cases, the intersectPolylines function may return the
desired result in a faster way.
PTS = polylineSelfIntersections(POLY);
Returns the position of self intersections of the given polyline.
PTS = polylineSelfIntersections(POLY, CLOSED);
Adds an options to specify if the polyline is closed (i.e., is a
polygon), or open (the default). CLOSED can be a boolean, or one of
'closed' or 'open'.
[PTS, POS1, POS2] = polylineSelfIntersections(POLY);
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).
[...] = polylineSelfIntersections(POLY, 'tolerance', TOL)
Specifies an additional parameter to decide whether two intersection
points should be considered the same, based on their Euclidean
distance.
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 =
Empty matrix: 0-by-2
polylineSelfIntersections(poly, 'closed')
ans =
10 10
See also
polygons2d, intersectPolylines, polygonSelfIntersections0001 function varargout = polylineSelfIntersections(poly, varargin) 0002 %POLYLINESELFINTERSECTIONS Find self-intersection points of a polyline. 0003 % 0004 % Computes self-intersections of a polyline, eventually specifying if 0005 % polyline is closed or open, and eventually returning position of 0006 % intersection points on polyline. 0007 % For common use cases, the intersectPolylines function may return the 0008 % desired result in a faster way. 0009 % 0010 % 0011 % PTS = polylineSelfIntersections(POLY); 0012 % Returns the position of self intersections of the given polyline. 0013 % 0014 % PTS = polylineSelfIntersections(POLY, CLOSED); 0015 % Adds an options to specify if the polyline is closed (i.e., is a 0016 % polygon), or open (the default). CLOSED can be a boolean, or one of 0017 % 'closed' or 'open'. 0018 % 0019 % [PTS, POS1, POS2] = polylineSelfIntersections(POLY); 0020 % Also return the 2 positions of each intersection point (the position 0021 % when meeting point for first time, then position when meeting point 0022 % for the second time). 0023 % 0024 % [...] = polylineSelfIntersections(POLY, 'tolerance', TOL) 0025 % Specifies an additional parameter to decide whether two intersection 0026 % points should be considered the same, based on their Euclidean 0027 % distance. 0028 % 0029 % 0030 % Example 0031 % % use a gamma-shaped polyline 0032 % poly = [0 0;0 10;20 10;20 20;10 20;10 0]; 0033 % polylineSelfIntersections(poly) 0034 % ans = 0035 % 10 10 0036 % 0037 % % use a 'S'-shaped polyline 0038 % poly = [10 0;0 0;0 10;20 10;20 20;10 20]; 0039 % polylineSelfIntersections(poly) 0040 % ans = 0041 % Empty matrix: 0-by-2 0042 % polylineSelfIntersections(poly, 'closed') 0043 % ans = 0044 % 10 10 0045 % 0046 % See also 0047 % polygons2d, intersectPolylines, polygonSelfIntersections 0048 % 0049 0050 % ------ 0051 % Author: David Legland 0052 % e-mail: david.legland@inra.fr 0053 % Created: 2009-06-15, using Matlab 7.7.0.471 (R2008b) 0054 % Copyright 2009 INRA - Cepia Software Platform. 0055 0056 % HISTORY 0057 % 2009/06/18 check bounding boxes before computing intersections 0058 0059 0060 %% Initialisations 0061 0062 % flag indicating whether the polyline is closed (polygon) or not 0063 closed = false; 0064 0065 % the tolerance for comparing positions based on distances 0066 tol = 1e-14; 0067 0068 % determine whether the polyline is open or closed 0069 if ~isempty(varargin) 0070 closed = varargin{1}; 0071 if ischar(closed) 0072 if strcmp(closed, 'closed') 0073 closed = true; 0074 varargin(1) = []; 0075 elseif strcmp(closed, 'open') 0076 closed = false; 0077 varargin(1) = []; 0078 end 0079 end 0080 end 0081 0082 % parse optional arguments 0083 while length(varargin) > 1 0084 pname = varargin{1}; 0085 if ~ischar(pname) 0086 error('Expect optional arguments as name-value pairs'); 0087 end 0088 0089 if strcmpi(pname, 'tolerance') 0090 tol = varargin{2}; 0091 else 0092 error(['Unknown parameter name: ' pname]); 0093 end 0094 varargin(1:2) = []; 0095 end 0096 0097 % if polyline is closed, ensure the last point equals the first one 0098 if closed 0099 if sum(abs(poly(end, :) - poly(1,:)) < tol) ~= 2 0100 poly = [poly; poly(1,:)]; 0101 end 0102 end 0103 0104 % arrays for storing results 0105 points = zeros(0, 2); 0106 pos1 = zeros(0, 1); 0107 pos2 = zeros(0, 1); 0108 0109 % number of edges 0110 nEdges = size(poly, 1) - 1; 0111 0112 0113 %% Main processing 0114 0115 % index of current intersection 0116 ip = 0; 0117 0118 % iterate over each couple of edge ( (N-1)*(N-2)/2 iterations) 0119 for iEdge1 = 1:nEdges-1 0120 % create first edge 0121 edge1 = [poly(iEdge1, :) poly(iEdge1+1, :)]; 0122 for iEdge2 = iEdge1+2:nEdges 0123 % create second edge 0124 edge2 = [poly(iEdge2, :) poly(iEdge2+1, :)]; 0125 0126 % check conditions on bounding boxes, to avoid computing the 0127 % intersections 0128 if min(edge1([1 3])) > max(edge2([1 3])) 0129 continue; 0130 end 0131 if max(edge1([1 3])) < min(edge2([1 3])) 0132 continue; 0133 end 0134 if min(edge1([2 4])) > max(edge2([2 4])) 0135 continue; 0136 end 0137 if max(edge1([2 4])) < min(edge2([2 4])) 0138 continue; 0139 end 0140 0141 % compute intersection point 0142 inter = intersectEdges(edge1, edge2, tol); 0143 0144 if sum(isfinite(inter)) == 2 0145 % add point to the list 0146 ip = ip + 1; 0147 points(ip, :) = inter; 0148 0149 % also compute positions on the polyline 0150 pos1(ip, 1) = iEdge1 - 1 + edgePosition(inter, edge1); 0151 pos2(ip, 1) = iEdge2 - 1 + edgePosition(inter, edge2); 0152 end 0153 end 0154 end 0155 0156 0157 %% Post-processing 0158 0159 % if polyline is closed, the first vertex was found as an intersection, so 0160 % we need to remove it 0161 if closed 0162 % identify the intersection between first and last edges using position 0163 % indices (pos1 < pos2 by construction) 0164 ind = pos1 == 0 & pos2 == size(poly,1)-1; 0165 points(ind,:) = []; 0166 pos1(ind) = []; 0167 pos2(ind) = []; 0168 end 0169 0170 % remove multiple intersections 0171 [points, I, J] = unique(points, 'rows', 'first'); %#ok<ASGLU> 0172 pos1 = pos1(I); 0173 pos2 = pos2(I); 0174 0175 % remove multiple intersections, using tolerance on distance 0176 iInter = 0; 0177 while iInter < size(points, 1) - 1 0178 iInter = iInter + 1; 0179 % for iInter = 1:size(points, 1)-1 0180 % determine distance between current point and remaining points 0181 inds = iInter+1:size(points, 1); 0182 dists = distancePoints(points(iInter,:), points(inds, :)); 0183 0184 % identify index of other points located in a close neighborhood 0185 inds = inds(dists < tol); 0186 0187 % remove redundant intersection points 0188 if ~isempty(inds) 0189 points(inds, :) = []; 0190 pos1(inds) = []; 0191 pos2(inds) = []; 0192 end 0193 end 0194 0195 0196 %% process output arguments 0197 0198 if nargout <= 1 0199 varargout{1} = points; 0200 0201 elseif nargout == 3 0202 varargout{1} = points; 0203 varargout{2} = pos1; 0204 varargout{3} = pos2; 0205 end