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@inrae.fr 0053 % Created: 2009-06-15, using Matlab 7.7.0.471 (R2008b) 0054 % Copyright 2009-2024 INRA - Cepia Software Platform 0055 0056 %% Initialisations 0057 0058 % flag indicating whether the polyline is closed (polygon) or not 0059 closed = false; 0060 0061 % the tolerance for comparing positions based on distances 0062 tol = 1e-14; 0063 0064 % determine whether the polyline is open or closed 0065 if ~isempty(varargin) 0066 closed = varargin{1}; 0067 if ischar(closed) 0068 if strcmp(closed, 'closed') 0069 closed = true; 0070 varargin(1) = []; 0071 elseif strcmp(closed, 'open') 0072 closed = false; 0073 varargin(1) = []; 0074 end 0075 end 0076 end 0077 0078 % parse optional arguments 0079 while length(varargin) > 1 0080 pname = varargin{1}; 0081 if ~ischar(pname) 0082 error('Expect optional arguments as name-value pairs'); 0083 end 0084 0085 if strcmpi(pname, 'tolerance') 0086 tol = varargin{2}; 0087 else 0088 error(['Unknown parameter name: ' pname]); 0089 end 0090 varargin(1:2) = []; 0091 end 0092 0093 % if polyline is closed, ensure the last point equals the first one 0094 if closed 0095 if sum(abs(poly(end, :) - poly(1,:)) < tol) ~= 2 0096 poly = [poly; poly(1,:)]; 0097 end 0098 end 0099 0100 % arrays for storing results 0101 points = zeros(0, 2); 0102 pos1 = zeros(0, 1); 0103 pos2 = zeros(0, 1); 0104 0105 % number of edges 0106 nEdges = size(poly, 1) - 1; 0107 0108 0109 %% Main processing 0110 0111 % index of current intersection 0112 ip = 0; 0113 0114 % iterate over each couple of edge ( (N-1)*(N-2)/2 iterations) 0115 for iEdge1 = 1:nEdges-1 0116 % create first edge 0117 edge1 = [poly(iEdge1, :) poly(iEdge1+1, :)]; 0118 for iEdge2 = iEdge1+2:nEdges 0119 % create second edge 0120 edge2 = [poly(iEdge2, :) poly(iEdge2+1, :)]; 0121 0122 % check conditions on bounding boxes, to avoid computing the 0123 % intersections 0124 if min(edge1([1 3])) > max(edge2([1 3])) 0125 continue; 0126 end 0127 if max(edge1([1 3])) < min(edge2([1 3])) 0128 continue; 0129 end 0130 if min(edge1([2 4])) > max(edge2([2 4])) 0131 continue; 0132 end 0133 if max(edge1([2 4])) < min(edge2([2 4])) 0134 continue; 0135 end 0136 0137 % compute intersection point 0138 inter = intersectEdges(edge1, edge2, tol); 0139 0140 if sum(isfinite(inter)) == 2 0141 % add point to the list 0142 ip = ip + 1; 0143 points(ip, :) = inter; 0144 0145 % also compute positions on the polyline 0146 pos1(ip, 1) = iEdge1 - 1 + edgePosition(inter, edge1); 0147 pos2(ip, 1) = iEdge2 - 1 + edgePosition(inter, edge2); 0148 end 0149 end 0150 end 0151 0152 0153 %% Post-processing 0154 0155 % if polyline is closed, the first vertex was found as an intersection, so 0156 % we need to remove it 0157 if closed 0158 % identify the intersection between first and last edges using position 0159 % indices (pos1 < pos2 by construction) 0160 ind = pos1 == 0 & pos2 == size(poly,1)-1; 0161 points(ind,:) = []; 0162 pos1(ind) = []; 0163 pos2(ind) = []; 0164 end 0165 0166 % remove multiple intersections 0167 [points, I, J] = unique(points, 'rows', 'first'); %#ok<ASGLU> 0168 pos1 = pos1(I); 0169 pos2 = pos2(I); 0170 0171 % remove multiple intersections, using tolerance on distance 0172 iInter = 0; 0173 while iInter < size(points, 1) - 1 0174 iInter = iInter + 1; 0175 % for iInter = 1:size(points, 1)-1 0176 % determine distance between current point and remaining points 0177 inds = iInter+1:size(points, 1); 0178 dists = distancePoints(points(iInter,:), points(inds, :)); 0179 0180 % identify index of other points located in a close neighborhood 0181 inds = inds(dists < tol); 0182 0183 % remove redundant intersection points 0184 if ~isempty(inds) 0185 points(inds, :) = []; 0186 pos1(inds) = []; 0187 pos2(inds) = []; 0188 end 0189 end 0190 0191 0192 %% process output arguments 0193 0194 if nargout <= 1 0195 varargout{1} = points; 0196 0197 elseif nargout == 3 0198 varargout{1} = points; 0199 varargout{2} = pos1; 0200 varargout{3} = pos2; 0201 end