Code covered by the BSD License

# geom2d

### David Legland (view profile)

13 Jun 2005 (Updated )

Geometry library for matlab. Performs geometric computations on points, lines, circles, polygons...

polylineSelfIntersections(poly, varargin)
```function varargout = polylineSelfIntersections(poly, varargin)
%POLYLINESELFINTERSECTIONS Find self-intersection points of a polyline
%
%   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).
%
%   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, polygonSelfIntersections
%

% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2009-06-15,    using Matlab 7.7.0.471 (R2008b)
% Copyright 2009 INRA - Cepia Software Platform.

%   HISTORY
%   2009/06/18 check bounding boxes before computing intersections

%% Initialisations

% determine whether the polyline is open or closed
closed = false;
if ~isempty(varargin)
closed = varargin{1};
if ischar(closed)
if strcmp(closed, 'closed')
closed = true;
elseif strcmp(closed, 'open')
closed = false;
end
end
end

% if polyline is closed, ensure the last point equals the first one
if closed
if sum(abs(poly(end, :) - poly(1,:)) < 1e-14) ~= 2
poly = [poly; poly(1,:)];
end
end

% arrays for storing results
points  = zeros(0, 2);
pos1    = zeros(0, 1);
pos2    = zeros(0, 1);

% number of vertices
Nv = size(poly, 1);

%% Main processing

% index of current intersection
ip = 0;

% iterate over each couple of edge ( (N-1)*(N-2)/2 iterations)
for i = 1:Nv-2
% create first edge
edge1 = [poly(i, :) poly(i+1, :)];
for j = i+2:Nv-1
% create second edge
edge2 = [poly(j, :) poly(j+1, :)];

% check conditions on bounding boxes, to avoid computing the
% intersections
if min(edge1([1 3])) > max(edge2([1 3]))
continue;
end
if max(edge1([1 3])) < min(edge2([1 3]))
continue;
end
if min(edge1([2 4])) > max(edge2([2 4]))
continue;
end
if max(edge1([2 4])) < min(edge2([2 4]))
continue;
end

% compute intersection point
inter = intersectEdges(edge1, edge2);

if sum(isfinite(inter)) == 2
% add point to the list
ip = ip + 1;
points(ip, :) = inter;

% also compute positions on the polyline
pos1(ip, 1) = i + edgePosition(inter, edge1) - 1;
pos2(ip, 1) = j + edgePosition(inter, edge2) - 1;
end
end
end

% if polyline is closed, the first vertex was found as an intersection, so
% we need to remove it
if closed
dist = distancePoints(points, poly(1,:));
[minDist, ind] = min(dist); %#ok<ASGLU>
points(ind,:) = [];
pos1(ind)   = [];
pos2(ind)   = [];
end

%% Post-processing

% process output arguments
if nargout <= 1
varargout{1} = points;

elseif nargout == 3
varargout{1} = points;
varargout{2} = pos1;
varargout{3} = pos2;
end
```

Contact us