Fast and Robust Curve Intersections

Computes intersection points of two curves.

Description

This function computes the (x,y) locations where two curves intersect. The curves can be broken with NaNs or have vertical segments. It is also very fast (at least on data that represents what I think is a typical application).

Acknowledgements

The author wishes to acknowledge the following in the creation of this submission:
Curve Intersect 2
This submission has inspired the following:
Fast and Robust Self-Intersections, Kirchhoff Vortex Contour Dynamics Simulation

• MATLAB release: MATLAB 7.0.4 (R14SP2)

[A. Silverstone]

Finally! A code that will work with vertical lines. Perfect timing on posting this too! Thanks!!!!

[M MA]

Well, the function meetpoint (http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=8997&objectType=FILE) already works with vertical lines... but the rate 2 is because I found one error at the first test!! using the example in meetpoint:
x1=[0 1.5 3 5 5];
y1=[0 1.5 1.5 2 1];
x2=[.2 .5 .5 1 1 2 2 5 5];
y2=[.2 .5 .8 1 1.3 1.3 1.8 1.8 .5];
[x,y,m1,m2]=meetpoint(x1,y1,x2,y2);
[x0,y0]=intersections(x1,y1,x2,y2);

The output of meetpoint is:
>> [x' y']
ans =
    0.5000 0.5000
    1.0000 1.0000
    1.3000 1.3000
    2.0000 1.5000
    4.2000 1.8000
    5.0000 1.8000

While the output of intersections is:
>> [x0 y0]
ans =
    0.5000 0.5000
    1.0000 1.0000
    1.3000 1.3000
    2.0000 1.5000
    4.2000 1.8000

So, it looks like its not handling vertical lines (both) correctly as meetpoint does!
Moreover, meetpoint also returns the slopes at the intersections!

[Douglas Schwarz]

M MA,
I'm sorry you felt you needed to give my function such a low rating, but it is not doing quite the same thing as yours (meetpoint). Comparing test suites is a little unfair under these circumstances. Nevertheless, I have attempted to make my function behave better when presented with pathological data. Speed and accuracy with real data are my priorities. I have just uploaded a new version -- it should be available soon.

[John D'Errico]

Doug has done a very nice job here, writing a code that is both fast and carefully constructed, watching out for degeneracies. Excellent explanations internally, in case you wish to know how he did it. Ignore the sour grapes review by M Ma. Or, if you prefer, consider the following comparison:

n = 100;
x1 = rand(1,n);
x2 = rand(1,n);
y1 = rand(1,n);
y2 = rand(1,n);

tic,[x0,y0]=intersections(x1,y1,x2,y2);toc
Elapsed time is 0.282314 seconds.

tic,[x0m,y0m]=meetpoint(x1,y1,x2,y2);toc
Elapsed time is 2.219815 seconds.

Both codes found the same points of intersections (2183 of them in total.) But intersections was considerably faster.

[Jos vdG]

Good job. Well coded and well documented.

[Peter Navé]

Excellent!

[Brian K]

Works very fine.

[Rickey Chen]

Excellent! Just what I am looking for. Thanks a lot!

[Chris Barichievy]

you're my hero.works like a dream; thanks very much

[Andy Dixon]

Does exactly what it says on the tin. Excellent!

[Ken Garrard]

A very nicely written and well documented function that is also efficient. I appreciate the error checking and handling of special cases. The comments make it straightforward to understand how the function works. I need to know which line segments contain intersection points and it was easy to add a couple of statements to return this information ( i(in_range) and j(in_range) ). I can confirm that the function is slower if A is changed from sparse to full. Does anyone know why?

This is exactly the sort of function that I hope to find whenever I browse the FEX.

[Eduardo GR]

The best intersection tool on FC, thanks.

[Gabriel Iniesta]

The function was working well, but i found no intersection for [X,Y]=intersections([1 6],[1 2],[6 8],[2 2]) when [6 2] actually is an intersection point. Other intersection tools do find this point....

[Douglas Schwarz]

Gabriel, thanks for the good example. The definition of intersection is debatable and your example only touches at the end points. Nevertheless, I have reverted to using \ in my algorithm instead of linsolve as it seems to be more accurate (at least in this case). This particular failure didn't occur in any of my test cases. The new version should be available shortly.

[Scott Miller]

This routine works well and quickly.

One general comment for users: "intersections" is a general curve intersection routine, not a polygon itersection routine, although it can be used for this as well. When testing polygons, add the first point of the boundary to the end of the boundary list - "intersections" does not add the last point automatically like "fill", nor is it intended to.

Scott

[Heikki Huttunen]

Works perfectly.

[Levi Kilcher]

This is a great function. Best intersection routine out there so far. No doubt.

One thing is missing though:
Returning the indices of the intersection. I think that Tim Meyers got this piece right with "intxy" ID:15239. He used floating point indices to indicate the index (read his description). However, his routine does not return all intersections (only the first), and I would guess it is not as fast/robust as this one.

I haven't looked through this function long enough to figure out how difficult such an addition would be. However, I will say that this routine would (IMHO) become MUCH more powerful/useful if some sort of index returning were included.

I'd be happy to work with you on implementing this if you gave me an idea where to start. I've done a bit of thinking about the "intersection problem", but clearly not as much as you.

Thanks,
Levi

[Douglas Schwarz]

Levi,

I'm not sure what use that would be, but it was trivial to implement and doesn't hurt the execution time so I have no objection to adding it. I have just uploaded a new version.

Doug

[Levi Kilcher]

Fantastic!

Thanks so much Doug.

FYI: the reason for wanting a index is for the case where you have a set of spatial vectors (x1,y1) and (x2,y2) but you also have a third set of data (i.e. t1) which is in one-to-one correspondence with (x1,y1). The indices which you have now included make it trivial to obtain the value of t1 where (x1,y1) and (x2,y2) intersect.

Hopefully this make sense.

Thanks again,
Levi

[W. H. Brave]

Excellent, thank you.

[Thomas]

Great code, works perfectly, and fast.

[Jeremy Bruch]

Doug,

Fanstastic job on this. The inclusion of indices is a great convenience as well. Thank you for sharing this code and continuing to improve on your work.

Jeremy

[Yuri Kotliarov]

Great function, works perfectly in almost all cases. However I found situation when it gives unpredictable result.

I draw 2 circles (using famous CIRCLE function from FEX), which touch each other.
clf
N=30;
hold on
H1=circle([0 0],5,N);
X1=get(H1,'XData');
Y1=get(H1,'YData');
H2=circle([0 10],5,N);
X2=get(H1,'XData');
Y2=get(H1,'YData');
[x,y]=intersections(X1,Y1,X2,Y2,0);
plot(x,y,'ro')
hold off
axis equal

It doesn't find intersection even with larger N and returns many warnings like
Warning: Matrix is singular to working precision.
or
Warning: Matrix is close to singular or badly scaled.
         Results may be inaccurate. RCOND = 1.559090e-016.

If I use N=10 or N=20 it actually finds intersection but in a weird place, specially at 20.

Thanks a lot anyway. 4 points.

[Douglas Schwarz]

Yuri,

If you zoom in on the supposed intersection point I think you'll find that there is no intersection at all. Remember, those circles are made up of line segments and so are actually polygons with N sides. For most values of N the two polygons won't even touch. In fact, the closest they can ever come is to touch at just one point and because of floating point round-off effects, my routine may or may not find that point.

Doug

[Yuri Kotliarov]

Douglas, thank you for comment.
That was exactly my point. At smaller N the polygons don't touch, but the function does find an intersection. If I use robust as true, much more intersection points on the 1st polygon are found.
At high N many points from two circles are actually very close, almost ideal touch. (Although I confirm there is a small gap at very close zoom in.) But at let's say N=1000, nothing found at robust=0, and 923 points all around the 1st circle at robust=1. Why is it?
Yuri

[Douglas Schwarz]

Yuri,

Please email me a mat-file and m-file that demonstrate what you are saying and I will investigate further.

Doug

[Douglas Schwarz]

It bothered me that I never resolved Yuri's issue so I just took another look at it and realized that his code has a bug. When he gets X2 and Y2 he is using the handle H1 instead of H2. After fixing that there is no problem.

[Yuri Kotliarov]

Doug, thanks and sorry for my stupid mistake and your time. I didn't test it yet, but believe that was the problem.

[Douglas Schwarz]

Yuri, no problem! I'm glad it's resolved and that you're still following this thread. -Doug

[Theresa]

Excellent file. It works perfectly, and extremely quickly. Just what I was looking for!

[Nik F]

 

[Olga Avila]

Excellent!
Found the link for this on the following thread <http://www.mathworks.de/matlabcentral/newsreader/view_thread/258260#672305&gt;
It was just what I needed :)