Code covered by the BSD License  

Highlights from
geom2d

image thumbnail

geom2d

by

 

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