Loren on the Art of MATLAB

Intersecting Lines 18

Posted by Loren Shure,

Today I am writing this post in collaboration with my friend and colleague, Lucio Cetto. Ever need to know if two line segments intersect? There are lots of ways to do this, and I thought I would show a few solutions here. I haven't put much thought into how to set up the data so I can vectorize the case for finding out which segments among many intersect and with which ones. At least for now :-)

Contents

Example

Let's start with with some simple line segments. The first column represents to x-values of the line segment, the second column represents the y-values.

line1 = [.5 1; 1 0];
line2 = [.5 0; 1 1];

Plot the data.

h = plot(line1(:,1),line1(:,2));
hold all
h(2) = plot(line2(:,1),line2(:,2));
set(h,'linewidth',2)
axis([0 1 0 1])

First Algorithm

Using equation y = mx+b, solve for x assuming 2 lines intersect. Then see if that x value is in the necessary range. Special cases: vertical lines (m==inf) and parallel lines (m1 == m2)

First find slopes and intercepts for both line segments. Here are slopes.

slope = @(line) (line(2,2) - line(1,2))/(line(2,1) - line(1,1));
m1 = slope(line1)
m2 = slope(line2)
m1 =
    -2
m2 =
     2

Here are the intercepts.

intercept = @(line,m) line(1,2) - m*line(1,1);
b1 = intercept(line1,m1)
b2 = intercept(line2,m2)
xintersect = (b2-b1)/(m1-m2)
yintersect = m1*xintersect + b1
b1 =
     2
b2 =
    -1
xintersect =
         0.75
yintersect =
          0.5

Plot results.

plot(xintersect,yintersect,'m*','markersize',8)
legend({'line1','line2','intersection'},'Location','West')
hold off

Now check that the intersection point is on the line segment and not past the end.

isPointInside = @(xint,myline) ...
    (xint >= myline(1,1) && xint <= myline(2,1)) || ...
    (xint >= myline(2,1) && xint <= myline(1,1));
inside = isPointInside(xintersect,line1) && ...
         isPointInside(xintersect,line2)
%
inside =
     1

Try with different lines now.

line1 = [.5 1; ...
        1 0];
line2 = [.5 0; ...
        0 -1];

Compute new intersection.

m1 = slope(line1);
m2 = slope(line2);
b1 = intercept(line1,m1);
b2 = intercept(line2,m2);
xintersect = (b2-b1)/(m1-m2)
yintersect = m1*xintersect + b1
xintersect =
         0.75
yintersect =
          0.5

Show the plot.

figure
h = plot(line1(:,1),line1(:,2));
hold all
h(2) = plot(line2(:,1),line2(:,2));
set(h,'linewidth',2)
legend(h,{'line1','line2'},'Location','West')
plot(xintersect,yintersect,'g*')
hold off
axis([-1 1 -1 1])

Is the intersect point on both line segments?

inside = isPointInside(xintersect,line1) && ...
         isPointInside(xintersect,line2)
inside =
     0

No. The intersect point lies on only one line segment, meaning that the two segments don't intersect.

There's still more to think about. What if the two lines segments have the same slope? There are 3 cases for this. First case: the segments are parallel and therefore no intersection. Second and third cases, they are part of the same line, in which case, we have to check to see if they overlap or not (in other words, they may intersect or not, for 2 more cases).

Here are two parallel lines with different values for their y-intercept.

line1 = [.5 1; ...
        1 0];
line2 = [.5 3; ...
        1 2];

Plot the line segments.

figure
h = plot(line1(:,1),line1(:,2));
hold all
h(2) = plot(line2(:,1),line2(:,2));
set(h,'linewidth',2)
legend(h,{'line1','line2'},'Location','West')
hold off

See that slopes are equal, and intercepts are not.

m1 = slope(line1);
m2 = slope(line2);
b1 = intercept(line1,m1);
b2 = intercept(line2,m2);
sameSlope = abs(m1-m2) < eps(m1)
differentIntercept = abs(b1-b2) > eps(b1)
isParallel = sameSlope && differentIntercept
sameSlope =
     1
differentIntercept =
     1
isParallel =
     1

Second Algorithm

Here's another algorithm for seeing if two lines intersect. The idea is to choose one line, and see if the end points from the other line lie on the same side. If they do, there's no way the lines have a point of intersection. If not, the second line might intersect the first one, or the point of intersection may fall outside the limits of the first line segment. A way to test that is to reverse the roles of lines 1 and 2 and do the test again.

The test for checking which side of the line a point falls on is simple. Plug in the values for the end points of line 2 into the equation for line 1. If they are both the same sign, then they both fall on the same side.

line1 = [.5 1; ...
        1 0];
line2 = [.5 0; ...
        0 -1];
figure
h = plot(line1(:,1),line1(:,2));
hold all
h(2) = plot(line2(:,1),line2(:,2));
set(h,'linewidth',2)
axis([0 1 -1 2])

Get slope and intercept of both lines.

m1 = slope(line1)
m2 = slope(line2)
b1 = intercept(line1,m1)
b2 = intercept(line2,m2)
m1 =
    -2
m2 =
     2
b1 =
     2
b2 =
    -1

Evaluate the x value of end-points of the second line using the equation y = m*x + b of the first line

lineEq = @(m,b, myline) m * myline(:,1) + b
yEst2 = lineEq(m1, b1, line2)
plot(line2(:,1),yEst2,'om')
lineEq = 
    @(m,b,myline)m*myline(:,1)+b
yEst2 =
     1
     2

Subtract the end values from the ends of line2 and we find that both values are positive, meaning we don't find that the segment straddles line2.

enderr = @(ends,line) ends - line(:,2)
errs1 = enderr(yEst2,line2)
enderr = 
    @(ends,line)ends-line(:,2)
errs1 =
     1
     3

Next reverse the roles of lines 1 and 2 and check to see if the end points of line 2 fall on the same side of line 1.

yEst1 = lineEq(m2, b2, line1)
plot(line1(:,1),yEst1,'or')
errs2 = enderr(yEst1,line1)
yEst1 =
     0
     1
errs2 =
    -1
     1

In this case we find that the differences between the predicted y values and the real y values of line 1 have different signs, which means that if line 2 was long enough there would exist an intersection.

To test for an actual intersection, we now just need to verify that the condition applies in both directions.

possibleIntersection =  sum(sign(errs1))==0 && sum(sign(errs2))==0
possibleIntersection =
     0

Ignored Issues and Other Algorithms

Here are a list of issues that this post does not solve:

  • Vertical line(s) (infinite slope)
        Can work around this by rotating the lines to avoid
        verticals (try random rotations, or perhaps rotate
        through a computed angle so largest angle is bisected
        vertically.
  • Degenerate line (single point)
  • Parallel lines with same intercept
  • More than two line segments

In addition, there is at least one more algorithm that we haven't talked about here, called Sweep Line Algorithm.

Ever Run Into This Problem Yourself?

How did you solve it? Let us know here.


Get the MATLAB code

Published with MATLAB® 7.12

18 CommentsOldest to Newest

I have run into an extension of this. I had two different tables each of which defined a set of line segments. I needed to find the intersection between the two sets.

I basically set up this exact problem with a layer of recursion. For each line segment from Table A, search for the intersect with each line segment from Table B. The recursion ended when the intersection fell within the domain for both line segments being tested.

This worked very well for a couple of years – literally hundreds of times – until I ran into a failure. The intersection actually fell on the domain boundary and the calculations in my code determined the intersection was outside the domain – by about 1 part in about 10^15. I had to cast the numbers to my desired resolution for the domain comparision to work in this case. It was a very painful lesson in the limitations of double-precision floating point!

I ran into a similar problem some time ago. In order to avoid the degenerate cases, I did use the second algorithm but I used the cross function to check on which side the endpoints of one segments were with respect to the other segment. This also opens the door to vectorization when you have more than two segments. Just build all the input for the cross function and perform all the cross products at once.

Dear Loren,
I think there’s an error on the eq. for the slope on the second algorithm. same eq. as for the first algorithm should be used.
the result, however, is correct, as line1=line1′ ans line2=line2′.

Hi Loren,

I scratched my idea in a paper and scanned it for you.
http://terpconnect.umd.edu/~cwchen/intersection%20problem.png

My idea is to parameterize the line segment. For example, PQ stands for a line segment. s is the parameter. If s=0, the point of interest is P. If s=1, the point of interest is Q. so s between [0,1] can stand for any point on the line segment.

Let’s do this parameterization for both line segments, and have s and t as the parameter for each segment, respectively. Then the question is

“Whether there exists a pair of [s,t] in the domain of [0,1] by [0,1]“

If the question on hand does, then we have intersection. Otherwise, we don’t.


If the intersection is not caused by two line “segments”, for example, just a line extending towards +/- infinity, my idea would still work. In detail,

first, pick two distinct points on the line first.
Second, parameterize the line as above with parameter s.
Third, replace s with tan(s) and domain of s is -pi/2 to pi/2

Then any point on the line can be mapped by assigning a s between -pi/2 to pi/2.

To test if two lines have intersection (of course your article has laid out the solution), define a domain [-pi/2,pi/2]^2 to see if there is a pair [s,t] locates in the domain.

To test if one line and one line segment has intersection,
define a domain [-pi/2,pi/2] by [0,1] to see if there is a pair [s,t] locates in the domain.

Translating and scaling [-pi/2,pi/2] to [0,1] is decorative. So I don’t want to go beyond.

I’ve read your blog for a while. I am excited I can help here. Hope my answer is satisfying.

Faithful reader

It looks like a lot of people were inspired by this post. Here is my contribution:
http://aprendtech.com/wordpress/?p=140

It uses complex variables to represent vectors in 2D space and describes segments by their endpoint vectors. With this, a point on a line is

R = P1 + s*D
where 
D = P2 - P1
and
0<=s<1

where all the capitalizedvariables are vectors represented as complex variables.

AFAIK it solves all the problems and issues you listed above. It can handle segments parallel to the y-axis with no fooling around, zero length segments, parallel segments, and more than one segment in the polyline.

I implemented fundamentally the same algorithm as Murat Erdem. This seems to me to be the natural way to parameterise the problem and avoid having to deal with special cases. The implementation in Matlab can be quite elegant:

function intersect = isintersect(line1, line2)
A = [line1(1,:) - line1(2,:); line2(2,:) - line2(1,:)]';
if det(A) == 0
    intersect = 0;
else
    mu = A \ (line2(2,:) - line1(2,:))';
    intersect = all(mu >= 0) && all(mu <= 1);
end

The problem could be rephrased as, does

  [A,-B;1,0;0,1]*s = [0;1;1]

have a solution with s≥0? (where A and B are the coordinate matrices of the vertices; for non-degenerate intersections we use s>0 instead). An advantage of this approach is that it works in higher dimensions and with convex polyhedrons. Besides linear-programming type solutions are there efficient general methods for this problem?

Said two line segments (A1-A2) and (B1-B2) where A1 A2 B1 B2 are the end points of the line segments, and A B are two points on the line segments respectively.

The question can also be thought of this way:

condition1= is Triangle A B1 B2 has zero area?
condition2= is Triangle B A1 A2 has zero area?

Is_intersection = condition1 && condition2

This idea also works for high dimension case.

I like Ben’s solution, but I have no clue to convert my own form to Ben’s.

If I am picking enough, to answer whether the intersection exists may be easier than to answer where the intersection is.

A=[0.5 1;1 0;0.5 0; 1 1];
% A = 4 (x,y) coordinates, order not important
in(4)=inpolygon(A(4,1),A(4,2),A([1 2 3],1),A([1 2 3],2));
in(3)=inpolygon(A(3,1),A(3,2),A([1 2 4],1),A([1 2 4],2));
in(2)=inpolygon(A(2,1),A(2,2),A([1 3 4],1),A([1 3 4],2));
in(1)=inpolygon(A(1,1),A(1,2),A([2 3 4],1),A([2 3 4],2));

is_intersection = all(in==0)

To use inpolygon, only 2D works.
To reduce 3D to 2D, one can test first whether the 4 points forming a tetrahedron have non-zero volume. If they do, they never have intersection. If they don’t, they are at least in a 2D plane. So one can then do test above.

The computation involved however may not be as small as other solutions. Didn’t do comparison yet. :)

If you wanted to be lazy, you could just call Java (this handles all of the special cases):

java.awt.geom.Line2D.linesIntersect(0.5, 1, 1, 0, 0.5, 0, 1, 1)

Similar to other algorithms described here, this function tests whether or not both endpoints of each line segment lie on opposite sides of the other line segment. What’s neat about the implementation is that it is able to perform the tests without division! It uses the signed area of each of the 4 triangles defined by the four points to determine the result. For example, if S is one line segment and (x1, y1) and (x2, y2) are the coordinates of the other line segment, then compute the signed area of S + (x1, y1) and the signed area of S + (x2, y2). If the signs of the areas are the same, then the two points lie on the same side of S, otherwise they are on opposite sides. There are some additional checks to run if one of the triangles has an area of zero, but that’s the basic idea.

Although I haven’t run any tests, I wouldn’t be surprised to learn this approach is quite fast, given the lack of floating point division.

You can view an implementation here:
http://developer.classpath.org/doc/java/awt/geom/Line2D-source.html

These are all very interesting and helpful.

I would be interested to see a similar post (and the comments) for finding the intersection between 2 lines in 3D space. Certainly not a straight-forward problem. Perhaps also the problem of finding a solution with some tolerance, or rather their closest point of intersection.

Thanks,
Val

These postings are the author's and don't necessarily represent the opinions of MathWorks.