How do you detect where two line segments intersect? [closed]

How do I determine whether or not two lines intersect, and if they do, at what x,y point?


Solution 1:

There’s a nice approach to this problem that uses vector cross products. Define the 2-dimensional vector cross product v × w to be vx wy − vy wx.

Suppose the two line segments run from p to p + r and from q to q + s. Then any point on the first line is representable as p + t r (for a scalar parameter t) and any point on the second line as q + u s (for a scalar parameter u).

Two line segments intersecting

The two lines intersect if we can find t and u such that:

p + t r = q + u s

Formulae for the point of intersection

Cross both sides with s, getting

(p + t r) × s = (q + u s) × s

And since s × s = 0, this means

t (r × s) = (qp) × s

And therefore, solving for t:

t = (qp) × s / (r × s)

In the same way, we can solve for u:

(p + t r) × r = (q + u s) × r

u (s × r) = (pq) × r

u = (pq) × r / (s × r)

To reduce the number of computation steps, it's convenient to rewrite this as follows (remembering that s × r = − r × s):

u = (qp) × r / (r × s)

Now there are four cases:

  1. If r × s = 0 and (q − p) × r = 0, then the two lines are collinear.

    In this case, express the endpoints of the second segment (q and q + s) in terms of the equation of the first line segment (p + t r):

    t0 = (qp) · r / (r · r)

    t1 = (q + sp) · r / (r · r) = t0 + s · r / (r · r)

    If the interval between t0 and t1 intersects the interval [0, 1] then the line segments are collinear and overlapping; otherwise they are collinear and disjoint.

    Note that if s and r point in opposite directions, then s · r < 0 and so the interval to be checked is [t1, t0] rather than [t0, t1].

  2. If r × s = 0 and (q − p) × r ≠ 0, then the two lines are parallel and non-intersecting.

  3. If r × s ≠ 0 and 0 ≤ t ≤ 1 and 0 ≤ u ≤ 1, the two line segments meet at the point p + t r = q + u s.

  4. Otherwise, the two line segments are not parallel but do not intersect.

Credit: this method is the 2-dimensional specialization of the 3D line intersection algorithm from the article "Intersection of two lines in three-space" by Ronald Goldman, published in Graphics Gems, page 304. In three dimensions, the usual case is that the lines are skew (neither parallel nor intersecting) in which case the method gives the points of closest approach of the two lines.

Solution 2:

FWIW, the following function (in C) both detects line intersections and determines the intersection point. It is based on an algorithm in Andre LeMothe's "Tricks of the Windows Game Programming Gurus". It's not dissimilar to some of the algorithm's in other answers (e.g. Gareth's). LeMothe then uses Cramer's Rule (don't ask me) to solve the equations themselves.

I can attest that it works in my feeble asteroids clone, and seems to deal correctly with the edge cases described in other answers by Elemental, Dan and Wodzu. It's also probably faster than the code posted by KingNestor because it's all multiplication and division, no square roots!

I guess there's some potential for divide by zero in there, though it hasn't been an issue in my case. Easy enough to modify to avoid the crash anyway.

// Returns 1 if the lines intersect, otherwise 0. In addition, if the lines 
// intersect the intersection point may be stored in the floats i_x and i_y.
char get_line_intersection(float p0_x, float p0_y, float p1_x, float p1_y, 
    float p2_x, float p2_y, float p3_x, float p3_y, float *i_x, float *i_y)
{
    float s1_x, s1_y, s2_x, s2_y;
    s1_x = p1_x - p0_x;     s1_y = p1_y - p0_y;
    s2_x = p3_x - p2_x;     s2_y = p3_y - p2_y;

    float s, t;
    s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / (-s2_x * s1_y + s1_x * s2_y);
    t = ( s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / (-s2_x * s1_y + s1_x * s2_y);

    if (s >= 0 && s <= 1 && t >= 0 && t <= 1)
    {
        // Collision detected
        if (i_x != NULL)
            *i_x = p0_x + (t * s1_x);
        if (i_y != NULL)
            *i_y = p0_y + (t * s1_y);
        return 1;
    }

    return 0; // No collision
}

BTW, I must say that in LeMothe's book, though he apparently gets the algorithm right, the concrete example he shows plugs in the wrong numbers and does calculations wrong. For example:

(4 * (4 - 1) + 12 * (7 - 1)) / (17 * 4 + 12 * 10)

= 844/0.88

= 0.44

That confused me for hours. :(

Solution 3:

The problem reduces to this question: Do two lines from A to B and from C to D intersect? Then you can ask it four times (between the line and each of the four sides of the rectangle).

Here's the vector math for doing it. I'm assuming the line from A to B is the line in question and the line from C to D is one of the rectangle lines. My notation is that Ax is the "x-coordinate of A" and Cy is the "y-coordinate of C." And "*" means dot-product, so e.g. A*B = Ax*Bx + Ay*By.

E = B-A = ( Bx-Ax, By-Ay )
F = D-C = ( Dx-Cx, Dy-Cy ) 
P = ( -Ey, Ex )
h = ( (A-C) * P ) / ( F * P )

This h number is the key. If h is between 0 and 1, the lines intersect, otherwise they don't. If F*P is zero, of course you cannot make the calculation, but in this case the lines are parallel and therefore only intersect in the obvious cases.

The exact point of intersection is C + F*h.

More Fun:

If h is exactly 0 or 1 the lines touch at an end-point. You can consider this an "intersection" or not as you see fit.

Specifically, h is how much you have to multiply the length of the line in order to exactly touch the other line.

Therefore, If h<0, it means the rectangle line is "behind" the given line (with "direction" being "from A to B"), and if h>1 the rectangle line is "in front" of the given line.

Derivation:

A and C are vectors that point to the start of the line; E and F are the vectors from the ends of A and C that form the line.

For any two non-parallel lines in the plane, there must be exactly one pair of scalar g and h such that this equation holds:

A + E*g = C + F*h

Why? Because two non-parallel lines must intersect, which means you can scale both lines by some amount each and touch each other.

(At first this looks like a single equation with two unknowns! But it isn't when you consider that this is a 2D vector equation, which means this is really a pair of equations in x and y.)

We have to eliminate one of these variables. An easy way is to make the E term zero. To do that, take the dot-product of both sides of the equation using a vector that will dot to zero with E. That vector I called P above, and I did the obvious transformation of E.

You now have:

A*P = C*P + F*P*h
(A-C)*P = (F*P)*h
( (A-C)*P ) / (F*P) = h

Solution 4:

I have tried to implement the algorithm so elegantly described by Jason above; unfortunately while working though the mathematics in the debugging I found many cases for which it doesn't work.

For example consider the points A(10,10) B(20,20) C(10,1) D(1,10) gives h=.5 and yet it is clear by examination that these segments are no-where near each other.

Graphing this makes it clear that 0 < h < 1 criteria only indicates that the intercept point would lie on CD if it existed but tells one nothing of whether that point lies on AB. To ensure that there is a cross point you must do the symmetrical calculation for the variable g and the requirement for interception is: 0 < g < 1 AND 0 < h < 1

Solution 5:

Here's an improvement to Gavin's answer. marcp's solution is similar also, but neither postpone the division.

This actually turns out to be a practical application of Gareth Rees' answer as well, because the cross-product's equivalent in 2D is the perp-dot-product, which is what this code uses three of. Switching to 3D and using the cross-product, interpolating both s and t at the end, results in the two closest points between the lines in 3D. Anyway, the 2D solution:

int get_line_intersection(float p0_x, float p0_y, float p1_x, float p1_y, 
    float p2_x, float p2_y, float p3_x, float p3_y, float *i_x, float *i_y)
{
    float s02_x, s02_y, s10_x, s10_y, s32_x, s32_y, s_numer, t_numer, denom, t;
    s10_x = p1_x - p0_x;
    s10_y = p1_y - p0_y;
    s32_x = p3_x - p2_x;
    s32_y = p3_y - p2_y;

    denom = s10_x * s32_y - s32_x * s10_y;
    if (denom == 0)
        return 0; // Collinear
    bool denomPositive = denom > 0;

    s02_x = p0_x - p2_x;
    s02_y = p0_y - p2_y;
    s_numer = s10_x * s02_y - s10_y * s02_x;
    if ((s_numer < 0) == denomPositive)
        return 0; // No collision

    t_numer = s32_x * s02_y - s32_y * s02_x;
    if ((t_numer < 0) == denomPositive)
        return 0; // No collision

    if (((s_numer > denom) == denomPositive) || ((t_numer > denom) == denomPositive))
        return 0; // No collision
    // Collision detected
    t = t_numer / denom;
    if (i_x != NULL)
        *i_x = p0_x + (t * s10_x);
    if (i_y != NULL)
        *i_y = p0_y + (t * s10_y);

    return 1;
}

Basically it postpones the division until the last moment, and moves most of the tests until before certain calculations are done, thereby adding early-outs. Finally, it also avoids the division by zero case which occurs when the lines are parallel.

You also might want to consider using an epsilon test rather than comparison against zero. Lines that are extremely close to parallel can produce results that are slightly off. This is not a bug, it is a limitation with floating point math.