--- /dev/null
+## Copyright (C) 2003-2011 David Legland <david.legland@grignon.inra.fr>
+## Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
+## All rights reserved.
+##
+## Redistribution and use in source and binary forms, with or without
+## modification, are permitted provided that the following conditions are met:
+##
+## 1 Redistributions of source code must retain the above copyright notice,
+## this list of conditions and the following disclaimer.
+## 2 Redistributions in binary form must reproduce the above copyright
+## notice, this list of conditions and the following disclaimer in the
+## documentation and/or other materials provided with the distribution.
+##
+## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
+## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
+## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+##
+## The views and conclusions contained in the software and documentation are
+## those of the authors and should not be interpreted as representing official
+## policies, either expressed or implied, of the copyright holders.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{pts} = } polylineSelfIntersections (@var{poly})
+## Find self-intersections points of a polyline
+##
+## Return the position of self intersection points
+##
+## 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
+## @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 =
+## 10 10
+## @end example
+##
+## @seealso{polygons2d, intersectPolylines, polygonSelfIntersections}
+## @end deftypefn
+
+function varargout = polylineSelfIntersections(poly, varargin)
+
+ ## Initialisations
+
+ # determine whether the polyline is 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
+ 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
+
+endfunction