Segment Intersection — Failing on Unknown Edge Cases

Question

I was recently at a programming competition. One of the problems was to determine whether two line segments are parallel, intersect, or do not intersect.

For all of the test cases I ran, my solution worked. However, when I submitted it, it was rejected (but no other information was given). I suspect there is either some edge case I am completely missing... or perhaps the judges messed up.

Input is entered as follows:

ax ay bx by cx cy dx dy

where the first segment is given by $$(a_x, a_y), (b_x, b_y)$$ and the second by $$(c_x, c_y), (d_x, d_y)$$

Some other rules/assumptions:

• We may assume that no segment has length zero; that is, each segment has distinct endpoints.

• Segments may share endpoint(s) with the other segment. Endpoints are included in the segment and thus are able to intersect (this also means that a segment whose endpoint lies within the other segment will intersect)

• Lines which both intersect and are parallel (either partial or complete overlap) should be considered parallel.

• All coordinates are integers on the open interval (-1000, 1000).

Here's a screenshot of the problem. Note that I'm ignoring the multiple tests option because I highly doubt that was the issue.

import java.util.Scanner;

public class QP1505
{
    public static void main(String[] args)
    {
        Scanner in = new Scanner(System.in);

        int ax = in.nextInt();
        int ay = in.nextInt();
        int bx = in.nextInt();
        int by = in.nextInt();
        int cx = in.nextInt();
        int cy = in.nextInt();
        int dx = in.nextInt();
        int dy = in.nextInt();

        in.close();

        /*
         * Cross multiply to check equal slopes This takes care of
         * +/- infinity as well as zero
         */
        if ((by - ay) * (dx - cx) == (dy - cy) * (bx - ax))
        {
            System.out.println("PARALLEL");
        }
        else
        {
            double m1 = (double) (by - ay) / (bx - ax);
            double m2 = (double) (dy - cy) / (dx - cx);

            // Find point of intersection
            double x, y;
            if (Double.isInfinite(m1))
            {
                x = ax;
                y = m2 * (x - cx) + cy;
            }
            else if (Double.isInfinite(m2))
            {
                x = cx;
                y = m1 * (x - ax) + ay;
            }
            else
            {
                x = (m1 * ax - ay - m2 * cx + cy) / (m1 - m2);
                y = m1 * (x - ax) + ay;
            }

            // Check bounds
            if (((x >= ax && x <= bx) || (x <= ax && x >= bx))
                    && ((y >= ay && y <= by) || (y <= ay && y >= by))
                    && ((x >= cx && x <= dx) || (x <= cx && x >= dx))
                    && ((y >= cy && y <= dy) || (y <= cy && y >= dy)))
            {
                System.out.println("INTERSECT");
            }
            else
            {
                System.out.println("DO NOT INTERSECT");
            }
        }

    }
}

And some test cases:

• Vertical Horizontal, Intersect

  0 0 0 5 -1 2 1 2
  INTERSECT
  

• Vertical Horizontal, No Intersect

  0 0 0 5 1 2 4 2
  DO NOT INTERSECT
  

• Horizontal Horizontal, Parallel

  0 0 5 0 1 1 5 1
  PARALLEL
  

• Vertical Regular, Intersect

  0 0 0 5 1 0 -1 4
  INTERSECT
  

• Vertical Regular, No Intersect

  0 0 0 5 1 0 4 6
  DO NOT INTERSECT
  

• Horizontal Regular, Intersect

  4 4 7 4 4 3 6 5
  INTERSECT
  

• Horizontal Regular, No Intersect

  2 2 8 2 3 3 9 6
  do not intersect
  

• Regular Regular, Intersect

  2 1 4 3 1 2 5 1
  INTERSECT
  

• Regular Regular, Parallel

  1 2 5 6 3 2 5 4
  PARALLEL
  

• Regular Regular, No Intersect

  1 2 5 6 3 2 8 5
  DO NOT INTERSECT
  

Let me know if anything else needs to be clarified (or any other test cases).

Answer 1

Failure cases

I ran your program and was able to find an endless number of failures. My technique was to have one line touch the origin (0,0) and then generate all lines that pass directly through the origin by simply choosing mirrored coordinates. I was hoping to find roundoff errors and I did. Here are just a small sampling of failure cases, all of which print "Do not intersect" instead of "intersect":

0 0 0 5 7 29 -7 -29
0 0 0 5 11 -30 -11 30
0 0 0 5 11 -15 -11 15
0 0 0 5 11 25 -11 -25
0 0 0 5 13 -30 -13 30
0 0 0 5 13 -15 -13 15
0 0 0 5 13 27 -13 -27
0 0 0 5 14 29 -14 -29
0 0 0 5 19 21 -19 -21
0 0 0 5 21 23 -21 -23
0 0 0 5 21 27 -21 -27
0 0 0 5 22 -30 -22 30
0 0 0 5 22 -15 -22 15
0 0 0 5 22 25 -22 -25
0 0 0 5 23 -26 -23 26
0 0 0 5 23 -13 -23 13
0 0 0 5 23 27 -23 -27
0 0 0 5 25 -29 -25 29
0 0 0 5 25 7 -25 -7
0 0 0 5 25 14 -25 -14
0 0 0 5 25 28 -25 -28
0 0 0 5 26 -30 -26 30
0 0 0 5 26 -15 -26 15
0 0 0 5 26 27 -26 -27
0 0 0 5 27 29 -27 -29
0 0 0 5 28 29 -28 -29
0 0 0 5 29 15 -29 -15

Floating point?

My personal opinion about solving this problem correctly is that if your coordinates are constrained to be within -1000..1000, then you can do the whole problem using integers and not use floating point at all.

What you would need to do is create a fraction class that holds a numerator and denominator. The fraction class needs to be able to: multiply, divide, add, subtract, compare (i.e. everything you are doing with your doubles). Since the range of coordinates is so small, you shouldn't have a problem with overflow. By doing this you will be safe from roundoff errors.

Answer 2

As it was already pointed out, computing the slope of the line segments (as a double) is problematic because the slope can be "infinite" for vertical segments, and because of rounding errors.

My suggestion is to use a different algorithm which does not compute the slope. If the coordinates are integers then all intermediate values are integers as well and no rounding errors can occur. (I would choose the same algorithm for floating point coordinates because it does not need any special cases for vertical or nearly vertical segments.)

The idea is to describe the line segment from \$(a_x, a_y) \$ to \$ (b_x, b_y) \$ in its parametric form $$ (x, y) = (a_x, a_y) + u (b_x - a_x, b_y - a_y) \, \quad \text{where } 0 \le u \le 1 \,. $$ Then the intersection of two line segments is a solution of the linear equation system $$ (b_x - a_x) u - (d_x - c_x) v = c_x - a_x \\ (b_y - a_y) u - (d_y - c_y) v = c_y - a_y $$ with \$ 0 \le u \le 1 \$ and \$ 0 \le v \le 1 \$. In the general case, the solution is given by $$ u = \frac {\Delta_u} \Delta \, , \quad v = \frac {\Delta_v} \Delta $$ with the determinants $$ \Delta = \begin{vmatrix} b_x - a_x & d_x - c_x \\ b_y - a_y & d_y - c_y \end{vmatrix} \, , \quad \Delta_u = \begin{vmatrix} c_x - a_x & d_x - c_x \\ c_y - a_y & d_y - c_y \end{vmatrix} \, , \quad \Delta_v = \begin{vmatrix} c_x - a_x & b_x - a_x \\ c_y - a_y & b_y - a_y \end{vmatrix} $$ If \$ \Delta = 0\$ then the line segments are parallel. Otherwise there is a unique solution \$u, v \$ to the linear equation system, and it easy to check if the solutions fall in the range \$ [0, 1] \$ without actually performing the division.

This leads to the following simple function:

static void checkIntersection(int ax, int ay, int bx, int by, int cx, int cy, int dx, int dy) {
    int det = (bx - ax)*(dy - cy) - (by - ay)*(dx - cx);
    if (det != 0) {
        /*
         * Lines intersect. Check if intersection point is on both segments:
         */
        int detu = (cx - ax)*(dy - cy) - (cy - ay)*(dx - cx);
        int detv = (cx - ax)*(by - ay) - (cy - ay)*(bx - ax);
        if (det < 0) {
            // Normalise to det>0 to simplify the following check.
            det = -det;
            detu = -detu;
            detv = -detv;
        }
        if (detu >= 0 && detu <= det && detv >= 0 && detv <= det) {
            System.out.println("INTERSECT");
        } else {
            System.out.println("NO NOT INTERSECT");
        }
    } else {
        /*
         * Lines are parallel (or identical):
         */
        System.out.println("PARALLEL");
    }
}

For integer coordinates in the range \$ -1000 \ldots 1000 \$, all the calculations can be done without overflow. Generally, for coordinates in a range \$ -M \ldots M \$, the computed values are bounded by \$ 8M^2\$ in absolute value, so you can choose the required type accordingly.

Answer 3

Abstraction

Code that's well abstracted is easier to maintain, ...and to debug. I'd start by defining what a segment is - it's an immutable data structure with 2 points, where each point is defined by X and Y values (I'm sure java has something like a Point class built-in somewhere).

So I would take this chunk:

    Scanner in = new Scanner(System.in);

    int ax = in.nextInt();
    int ay = in.nextInt();
    int bx = in.nextInt();
    int by = in.nextInt();
    int cx = in.nextInt();
    int cy = in.nextInt();
    int dx = in.nextInt();
    int dy = in.nextInt();

    in.close();

And extract it into its own getSegments method, so that this:

if ((by - ay) * (dx - cx) == (dy - cy) * (bx - ax))

Can look something like this:

if (isParallel(segment1, segment2))

And then you would have isParallel on its own, and you could write a dozen unit tests just for that method, covering every possible edge case.

---

Then you can extract another findIntersectionPoint method here, which would return a Point:

    double m1 = (double) (by - ay) / (bx - ax);
    double m2 = (double) (dy - cy) / (dx - cx);

    // Find point of intersection
    double x, y;
    if (Double.isInfinite(m1))
    {
        x = ax;
        y = m2 * (x - cx) + cy;
    }
    else if (Double.isInfinite(m2))
    {
        x = cx;
        y = m1 * (x - ax) + ay;
    }
    else
    {
        x = (m1 * ax - ay - m2 * cx + cy) / (m1 - m2);
        y = m1 * (x - ax) + ay;
    }

And then you can cover that method with another dozen unit tests, and - assuming the "check bounds" part is also extracted into its own method, simplify the calling code to this:

Point intersection = findIntersectionPoint(segment1, segment2);
if (intersects(intersection, segment1, segment2))
{
    // intersect
}
else
{
    // no intersect
}

And then you can bombard findIntersectionPoint and intersects with another dozen unit tests to cover all edge cases.

---

The key point in this answer is granularity: you have one single method that seems to work, but it's doing many things that can't be tested individually - so even if you do end up figuring out a failing test, you wouldn't know exactly where the problem is and where to start looking.

Defining a Segment object might be over-the-top, but I find it simplifies the code, at least how it reads: seeing segment2.point1.x instead of cx makes it much easier to mentally map the values.. at least in my poor little fried brain.

---

One thing that strikes me, as a c# programmer, is that your curly braces are... c#-style. Typically java code would look something like this:

if (Double.isInfinite(m1)) {
    x = ax;
    y = m2 * (x - cx) + cy;
} else if (Double.isInfinite(m2)) {
    x = cx;
    y = m1 * (x - ax) + ay;
} else {
    x = (m1 * ax - ay - m2 * cx + cy) / (m1 - m2);
    y = m1 * (x - ax) + ay;
}

Answer 4

Okay I found some actual examples. with integer coordinates.

EXAMPLE #1 - Wraparound

Let's try the most extreme case, where one of the segments is: (INTEGER.MIN_VALUE,0) (Integer.MAX_VALUE,1). The slope should be positive and nearly zero, but calculation fails spectacularly. In an effort to avoid integer division, the code casts the numerator to a double, but the problem is that the subtractions are still being performed using int arithmetic:

        double m1 = (double) (by - ay) / (bx - ax);
        m1 = (double)(1-0)/(Integer.MAX_VALUE-Integer.MIN_VALUE);

Here, Java wraps around:

        m1 = (1.0)/(-1);
        m1 = -1.0;

I don't think we need to continue here; the calculation of x will obviously be wrong. The wraparound could be avoided by making sure that the individual elements are each cast to double or java.Math.BigInteger, so that all arithmetic is performed with doubles or BigIntegers, but there are still problems with that approach -- see below.

EXAMPLE #2 - Roundoff error

In this example we make sure that wraparound doesn't happen, either by choosing an interval that doesn't overflow, or by casting to double. There's still a problem. Say we have (0,0) (Integer.MAX_VALUE,1) and (0,0) (Integer.MAX_VALUE-1,1).

        double m1 = (double) (by - ay) / (bx - ax);
        m1 = (double)(1 - 0)/(Integer.MAX_VALUE-0);
        m1 = (double)(1)/Integer.MAX_VALUE;
        m1 = 1.0/Integer.MAX_VALUE;
        m1 = 0b0.0000000000000000000000000000001000000000000000000000000000000100000000000000000000000;
        double m2 = (double) (dy - cy) / (dx - cx);
        m2 = (double)(1 - 0)/(Integer.MAX_VALUE-1-0);
        m1 = (double)(1)/(Integer.MAX_VALUE-1);
        m2 = 1.0/(Integer.MAX_VALUE-1);
        m2 = 0b0.0000000000000000000000000000001000000000000000000000000000001000000000000000000000000;

These slopes are small, and very close to each other but not quite identical. They do differ in fact, but the difference is calculated as 0b0.00000000000000000000000000000000000000000000000000000000000001 = 2^-62 when really it's slightly larger than that (it's 2^-62 + 2^-92 + ...). So, m1-m2 is either slightly too large or too small, which causes x to be either slightly too large or too small. So when the point of intersection is the endpoint of one of the segments, the algorithm determines the lines don't intersect, when in fact they do.

You can verify the binary floating point arithmetic using http://www.exploringbinary.com/binary-calculator/ . Do the calculation once to find out where the most significant digit is, then add on 53 more bits.

I also notice the parallel check needs wraparound protection, both on the subtractions and on the multiplications; with a little work I could find two segments that report as parallel when they're not, or report as non parallel when they are. If you cast everything to double you still run the risk of the 53 bits of precision not being enough (think about multiplying two 30-bit quantities). I think you may need to use java.math.BigDecimal or java.math.BigInteger here.

Answer 5

You're using a dangerous method to obtain the point of intersection. Your method works fine when one of the slopes is calculated as Infinite but it's numerically unstable if the two slopes differ by many orders of magnitude (more than the 15-17 decimal digits that a Java double supports). This can't happen because all the coordinates are integers. Even the most extreme points can't cause roundoff error.

The problem is with the slopes. When you calculate m1 and m2 the subtractions and divisions may not be exact, which introduces some error into the slopes themselves. This is then compounded when the two slopes are then subtracted from each other m1 - m2 in the quotient where you are calculating x. It's quite possible that this subtraction will remove all the signal and be left with only noise, which will make the calculated x completely invalid. While I wish this was true, because all the coordinates are ints the subtractions are exact, and the divisions don't introduce that much error.

How about the first segment being (0,1), (10^{17},10^{17}) and the second segment being (1,0), (10^{17},10^{17}). The lines are not parallel; they intersect at x == 1e17, y == 1e17. Since your test for parallel is correct, it will not print "parallel" -- so far so good. Unfortunately, due to the limits of Java's double, both m1 and m2 will be calculated as 1.0000000000. Your program will crash when it tries to compute x due to the divide by zero.

That doesn't work, outside range for int plus they'd be detected as parallel. No bounty for me.