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, Fast Line Segment Intersection, Gear shift map for Automatic Transmission in S mode, polymorph

• MATLAB release: MATLAB 7.8 (R2009a)

[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>;
It was just what I needed :)

[Sören Sieberling]

Perfect

[Philips]

An excellent function. Met my needs 100%. Thanks for making life easier for me.

[Rodrigo Portugal]

A really really excellent function. I'm using all the time. However I have a suggestion. When the curves do not interecept, like in this example
[xi,yi,i,j] = intersections([0,1],[0,0,],[1,0],[1,1])
instead to produce an error, it would be better, in my opinion, if it produces empty matrices as an output. This is partially done when the function is called with two outputs, such as,
[xi,yi] = intersections([0,1],[0,0,],[1,0],[1,1])
but not with four. I give you five points anyway!
Rodrigo

[Douglas Schwarz]

Rodrigo, I agree. I just fixed it and uploaded a new version -- it should appear soon. Thanks for the report.

[Rodrigo Portugal]

Douglas,
there is still a little bug related to mismatching dimensions of empty matrices. When we try, for instance,
[xi,yi,i,j] = intersections([0,1],[0,0,],[1,0],[1,1])
that's OK. However, if we try
[xi,yi,i,j] = intersections([0,1],[0,0,],[1,0.5,0],[1,0.5,1])
there is the error "Matrix dimensions should agree", because empty matrices should have compatible dimensions, even though they are empty entities.

So I have look at the code, and I have inserted the following code
i = i(:); j = j(:);
just after "i" and "j" definitions at line 165. In this way it seems that this little annoying bug was removed. Could you check it please?

Personally (and strangely) I am very pleased about this bug, because reading the help on empty matrices I've learned that they have a strong operational structure underneath. As a matter of fact, they behave like non-empty matrices, when we try to make math operations on them!

Congratulations and thanks for the prompt response.

Rodrigo

[Douglas Schwarz]

Arrgh! I actually did exactly that and then removed it thinking (mistakenly) that it wasn't needed. I ended up using reshape(i,[],1) rather than i(:) because it is documented to force the second dimension to be 1 even if i starts out as 0-by-0. i(:) does the right thing today, but maybe not in the future. Thanks again, Rodrigo! (New version just uploaded.)

[Joshua Arnott]

 

[Shadab Khan]

Excellent! It does exactly what is needed and saved so much of time. I am using this function heavily, it saved 300 lines of my code which was inefficient in finding the intersection points.

A suggestion to all: Ignore the post by M MA, this function works perfect, and is fast as stated by Douglas.

Thanks! And i hope you will help us all with such great functions in the future.

Shadab

[Rui Xavier]

Hi douglas!! i am currently doing a final year project on ear recognition.the idea of finding the feature vector to be trained involves finding the distances of a centroid from the intersection points(lines have been drawn at different angles to intersect a edge detected ear image)..but we are unable to find the coordinates of the intersections.hoping u can help us out asap

thanking u in anticipation
Rui

[Douglas Schwarz]

Rui, please email me directly with details of your problem and I will see if I can help.
Doug

[Rui Xavier]

Ya thanks Doug. Didnt have your email address so contacted you with a private message on matlab central itself. If you can give me you email i can mail you our code so far & ear images as an attachment so you can get a better idea. Thanks

[U. Murat Erdem]

Excellent work!

I have submitted a code doing similar task but faster. You can find it here:

http://www.mathworks.com/matlabcentral/fileexchange/27205

I would appreciate your feedback.

[Truc Phan]

Thanks for sharing.

However, can you please tell us the major difference and/or advantage between this function and the polyxpoly function in Mapping Toolbox. Thanks

[Andrey]

Reliable and robust, but it takes 6.6 seconds for two 10,000 point-functions with ~20 intersections, with the robust flag off. May be I can figure out how to make it less robust, but a bit faster.

[Evgeny Pr]

2Douglas Schwarz
Thank you! Thank you! Thank you!

2Truc Phan
1. Not all have the Mapping Toolbox.
2. Polyxpoly does not work when you want to find the intersection curve itself (loops in curve).

[David]

May I suggest an improvement to your code.

Currently the code checks for intersecting domains and then for intersecting domain checks if the lines actually cross. This op takes 16 flops (that I count) per pair of lines, with then the LU decomp and back substitution for a 4 by 4 system of equations where the domains do cross.

However if the intersection condition is stated in terms of L1 = O1 + k1*D1 and L2 = O2 + k2*D2 where L1 == L2, you are left with an equation in 2 unknowns (k1 and k2). This involves ~7 operations to find k1 and k2, with four more to find x0 and y0.

N.B. Your implementation of this was far faster than my one using matlab classes (of the order of 3000-4000 times faster) :).

Cheers

David

[Rebecca]

The function seems great but I get an error about the try/catch loop syntax for my version of matlab.
changing line 218 to "catch, err" avoids the error.

[Manuel]

Thanks for the file. Saved me lines of code and hours of computing time.

[Diarmuid]

Excellent code ... well written, fast and above all it works!!! Has saved me a huge amount of coding time!

[Laura]

I have asked my question wrongly. The intersection is outside the interval of the selected points:
e.g;

C1=(627,338),(619,288 )
C2= (621,351), (627,339)

because this gives empty for me and it should intersect

[Laura]

Dear Douglas,

My question: is it possible to obtain with you code the intersection of 2 curves if their intersection is outside the endpoints given?
thank you!
Laura

[Douglas Schwarz]

Laura,

You have defined two curves consisting of just one line segment each. These two segments do not intersect so that is why my function returns empty vectors.

In order to find the intersection of the two infinite straight lines defined by the segments you can set

x1=[627 619];
y1=[338 288];
x2=[621 627];
y2=[351 339];

and then follow the calculations in the comments in intersections.m. Set up the A and B matrices and then solve with T = A\B. Since you are looking for an intersection outside the endpoints you can ignore t1 and t2 and just look at x0 and y0.

Doug

[NURUL]

Hi, I have used this code to find the coordinates intersection between 2 vectors, where X1, Y1 is my latitude and longitude for vector 1, and X2,Y2 is my latitude and longitude for vector 2.
Previously this function works very well in my data, but now, I am not so sure why it is sometimes cannot find the intersections, even it is implemented on the same datasets. I am really wondering why?

[Douglas Schwarz]

NURUL, please send me an example.

[Jaime]

Thanks a lot bruder !!!

[Vooi Leng Tan]

Thanks a million!

[Tom]

Thank you so much for sharing this wonderful code with us.

One small issue though, I found a case where instead of detecting 3 intersections it finds 5. Perhaps it is pathological but let me know if you would like more info.

[Daniel]

Using in combination with polymorph(), thanks!

[Jack]

Hi,

I've got a pollygon that encloses the origin and I have 2 points P1 and P2 outside the pollygon such that the line P1P2 intersects the pollygon in 2 points, I1 and I2, that I determine using the 'intersections' function. Then I rotate my line P1P2 by 10 degrees 36 times and I determine the intersection points with the pollygon for each angle. The results are stored in 2 vectors, x0 and y0, each with 2 components. Now,what I am trying to figure out here is how does the function choose to allocate the numerical values of the coordinates of the 2 points to x0(1) and x0(2) respectively?More explicitly, what does x0(1) represent, the x coordinate of which one of the 2 intersection points? Is it random? Cheers!

[Douglas Schwarz]

Jack,

No, it's not random.

For the NON-ROBUST option:
Number the line segments of your polygon from 1 to N, as defined by the order of the x- and y-coordinates. For each P1P2 line segment, the intersections are found in order of the polygon segments. In other words, if P1P2 crosses polygon segments 12 and 25, the intersection with segment 12 will be first followed by the intersection with segment 25, because 12 < 25.

For the ROBUST option:
Because unique() is run on the intersections, they end up being sorted in increasing order of the x-coordinate (with any tie broken by sorting y).

Doug

[Douglas Schwarz]

Jack,

It just occurred to me that you can use the 3rd and 4th output arguments to tell you in which line segments the intersections lie. See the help for more info.

Doug

[Brock]

Very nice work. The option to identify index pairs between which the intersection points lie is a nice touch. This saves me several hours of work, and the code has worked great on my test cases so far. Thanks so much!

[Floris]

Very nice work. Thanks for sharing.

I agree that the index identification is a very useful extra feature.

[KYAW KYAW]

Dear Sir/Madam,

I am looking for the intersection of dynamic waves and base line. SO far, I could manage to get few of the intersection points using ratio (Y data of dynamic wave/ Y data of base line) for whole of samples. Pls advise me how should I continue since I was stuck in here so long.

I used your function but mine is dynamic wave form and base straing line so the X -axis is not constant line. So the array is only (1900x1) for wave form and base line. The Y values of wave form is changing but the base line is seems like constant. Let say, there is 64 intersection points visuality but what i can get is only 47 but some are not accutare positions.

Thanks and best regards
Kyaw

[Harrison]

 

[Björn]

Really good, faster and as useful as polyxpoly.

[David]

Excellent work

[Paul Koppen]

This is an absolutely fantastic piece of code. 5* well deserved. Thank you!

That said I really regret writing this but...

To me the returned values for I and J feel awkward. Because X0 and Y0 return the intersection points in terms of the original space in which the curves lie (X,Y), I would expect the second set of returned values to give coordinates in terms of the curves themselves. So if a curve has (Euclidean) length 4.8, then I would expect 2.4 to be halfway the curve, not at 4 tenth of the second line segment.

This way we could also compute derivatives and integrals along a curve and interpolate these values simply by looking op the value at I (or J for that matter). Furthermore, the distance along the curve is a nice linear function, whereas the currently returned values for I and J are piecewise linear and non-differentiable it the points.

I understand that such a change of I and J would not be backwards compatible, and also it might be slightly (but hardly) more involved to compute, so I see the big arguments against. Nevertheless, I would be very happy if you would at least give it a thought :-)

[Douglas Schwarz]

Paul,

Thank you for your comments. First of all, you should know that the definitions for I and J were specifically requested by users -- they were not included at all in early versions of the program. I have no opinion about them since I don't use them myself. Certainly, backwards compatibility is important so they are not going to change.

I would encourage you to write a simple wrapper function that runs intersections and then converts its outputs into whatever you require. I hesitate to do that for you since I'm not completely sure what it is you need, but feel free to contact me via email if you need help with it.

Doug

[Murat Shagirov]