This function computes the locations where a curve self-intersects in a fast and robust way. The curve can be broken with NaNs or have vertical segments. Segments of the curve involved in each of the self-interesections are also provided.

Example of use:
N=201;
th=linspace(-3*pi,4*pi,N);
R=1;
x=R*cos(th)+linspace(0,6,N);
y=R*sin(th)+linspace(0,1,N);
t0=clock;
[x0,y0,segments]=selfintersect(x,y)
etime(clock,t0)
plot(x,y,'b',x0,y0,'.r');
axis ('equal'); grid