Line Segments Intersection(intersection Point)

Question

I have created a function to calculate the intersection point of two line segment .

Unfortunantly the code below dosen't work if one of the segment is verticale

    public static Point intersection(Segment s1, Segment s2) {
    double x1 = s1.getP1().getX();
    double y1 = s1.getP1().getY() ;
    double x2 = s1.getP2().getX();
    double y2 = s1.getP2().getY() ;
    double x3 = s2.getP1().getX();
    double y3 = s2.getP1().getY();
    double x4 = s2.getP2().getX();
    double y4 = s2.getP2().getY();

    double d = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4);
    if (d == 0) {
        return null;
    }
    double xi = ((x3 - x4) * (x1 * y2 - y1 * x2) - (x1 - x2) * (x3 * y4 - y3 * x4)) / d;
    double yi = ((y3 - y4) * (x1 * y2 - y1 * x2) - (y1 - y2) * (x3 * y4 - y3 * x4)) / d;
    Point p = new Point(xi, yi);
    if (xi < Math.min(x1, x2) || xi > Math.max(x1, x2)) {
        return null;
    }
    if (xi < Math.min(x3, x4) || xi > Math.max(x3, x4)) {
        return null;
    }
    return p;
}

the problem when i have a vertical line segment , this formula

double d = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4);

is equal to 0 and the method return null.

How can I handle this exception.

Thank you

Answer 1

// calculate the intersection of 2 line segments
// segment 1 = points A and B
// segment 2 = points C and D
// if there is an intersection return the point in pt
// from https://en.wikipedia.org/wiki/Line%E2%80%93line_intersection
//
//      /      \     /         \          /      \     /         \
//      |  x1  |     | x2 - x1 |          |  x3  |     | x4 - x3 |
// L1 = |  --  | + t | ------- |,    L2 = |  --  | + u | ------- |
//      |  y1  |     | y2 - y1 |          |  y3  |     | y4 - y3 |
//      \      /     \         /          \      /     \         /
//
//
//      |  x1 - x3   x3 - x4  |
//      |  y1 - y3   y3 - y4  |     (x1 - x3)(y3 - y4) - (y1 - y3)(x3 - x4)
// t  = -----------------------  =  ---------------------------------------
//      |  x1 - x2   x3 - x4  |     (x1 - x2)(y3 - y4) - (y1 - y2)(x3 - x4)
//      |  y1 - y2   y3 - y4  |
//
//
//      |  x1 - x3   x1 - x2  |
//      |  y1 - y3   y1 - y2  |     (x1 - x3)(y1 - y2) - (y1 - y3)(x1 - x2)
// u  = -----------------------  =  ---------------------------------------
//      |  x1 - x2   x3 - x4  |     (x1 - x2)(y3 - y4) - (y1 - y2)(x3 - x4)
//      |  y1 - y2   y3 - y4  |
//
// (Px,Py) = (x1 + t(x2 - x1), y1 + t(y2 - y1))
//
// NOTE:
//    denominators are the same for t and u and must not be 0
//    t and u must be in the range 0 to 1
//
bool calc_intersection(const point& A, const point& B, const point& C, const point& D, point& pt)
{
   const double compare_tolerance = .00001;

    pt = { 0,0 };

    const auto x1 = A.x;
    const auto y1 = A.y;
    const auto x2 = B.x;
    const auto y2 = B.y;
    const auto x3 = C.x;
    const auto y3 = C.y;
    const auto x4 = D.x;
    const auto y4 = D.y;

    // this is done for simplifying terms
    const auto x1_x2 = x1 - x2;
    const auto x1_x3 = x1 - x3;
    const auto x2_x1 = x2 - x1;
    const auto x3_x4 = x3 - x4;

    const auto y1_y2 = y1 - y2;
    const auto y1_y3 = y1 - y3;
    const auto y2_y1 = y2 - y1;
    const auto y3_y4 = y3 - y4;

    const auto denominator = x1_x2 * y3_y4 - y1_y2 * x3_x4;
    if (abs(denominator) < compare_tolerance)
        return false;

    const auto t = (x1_x3 * y3_y4 - y1_y3 * x3_x4) / denominator;
    if (t < 0 || t > 1)
        return false;

    const auto u = (x1_x3 * y1_y2 - y1_y3 * x1_x2) / denominator;
    if (u < 0 || u > 1)
        return false;

    pt = point(x1 + t * x2_x1, y1 + t * y2_y1);
    return true;
}

Answer 2

Here's my answer. I have tested it to be accurate by making a loop that checks that the answer it gives is the same as the one given by the Boost geometry library and they agree on each test, though the one that I have written below is much much faster than the one in Boost. The test makes every possible line segment where x is and integer in [-3,2] and y is an integer in [-3,2], for all possible pairs of line segments.

The code below considers line segments that touch at endpoints to be intersecting. T-shaped intersections are also considered intersecting. The code is in c++ but would be easily adaptable to any language. It's based on a different stackoverflow answer but that answer did not handle endpoints correctly.

It uses the cross product method, which can report if a point is to the left or to the right of a given ray.

There are some optimizations to be made in the math but doing them showed no performance improvement over compilation with g++ -O2 and sometime the performance even decreased! The compiler is able to do those optimizations so I prefered to leave the code readable.

// is_left(): tests if a point is Left|On|Right of an infinite line.
//    Input:  three points p0, p1, and p2
//    Return: >0 for p2 left of the line through p0 and p1
//            =0 for p2 on the line
//            <0 for p2 right of the line
//    See: Algorithm 1 "Area of Triangles and Polygons"
//    This is p0p1 cross p0p2.
extern inline coordinate_type_fp is_left(point_type_fp p0, point_type_fp p1, point_type_fp p2) {
  return ((p1.x() - p0.x()) * (p2.y() - p0.y()) -
          (p2.x() - p0.x()) * (p1.y() - p0.y()));
}

// Is x between a and b, where a can be lesser or greater than b.  If
// x == a or x == b, also returns true. */
extern inline coordinate_type_fp is_between(coordinate_type_fp a,
                                            coordinate_type_fp x,
                                            coordinate_type_fp b) {
  return x == a || x == b || (a-x>0) == (x-b>0);
}

// https://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect
extern inline bool is_intersecting(const point_type_fp& p0, const point_type_fp& p1,
                                   const point_type_fp& p2, const point_type_fp& p3) {
  const coordinate_type_fp left012 = is_left(p0, p1, p2);
  const coordinate_type_fp left013 = is_left(p0, p1, p3);
  const coordinate_type_fp left230 = is_left(p2, p3, p0);
  const coordinate_type_fp left231 = is_left(p2, p3, p1);

  if (p0 != p1) {
    if (left012 == 0) {
      if (is_between(p0.x(), p2.x(), p1.x()) &&
          is_between(p0.y(), p2.y(), p1.y())) {
        return true; // p2 is on the line p0 to p1
      }
    }
    if (left013 == 0) {
      if (is_between(p0.x(), p3.x(), p1.x()) &&
          is_between(p0.y(), p3.y(), p1.y())) {
        return true; // p3 is on the line p0 to p1
      }
    }
  }
  if (p2 != p3) {
    if (left230 == 0) {
      if (is_between(p2.x(), p0.x(), p3.x()) &&
          is_between(p2.y(), p0.y(), p3.y())) {
        return true; // p0 is on the line p2 to p3
      }
    }
    if (left231 == 0) {
      if (is_between(p2.x(), p1.x(), p3.x()) &&
          is_between(p2.y(), p1.y(), p3.y())) {
        return true; // p1 is on the line p2 to p3
      }
    }
  }
  if ((left012 > 0) == (left013 > 0) ||
      (left230 > 0) == (left231 > 0)) {
    if (p1 == p2) {
      return true;
    }
    return false;
  } else {
    return true;
  }
}

The test code:

BOOST_AUTO_TEST_CASE(small, *boost::unit_test::disabled()) {
  for (double x0 = -3; x0 < 3; x0++) {
    for (double y0 = -3; y0 < 3; y0++) {
      for (double x1 = -3; x1 < 3; x1++) {
        for (double y1 = -3; y1 < 3; y1++) {
          for (double x2 = -3; x2 < 3; x2++) {
            for (double y2 = -3; y2 < 3; y2++) {
              for (double x3 = -3; x3 < 3; x3++) {
                for (double y3 = -3; y3 < 3; y3++) {
                  point_type_fp p0{x0, y0};
                  point_type_fp p1{x1, y1};
                  point_type_fp p2{x2, y2};
                  point_type_fp p3{x3, y3};
                  linestring_type_fp ls0{p0,p1};
                  linestring_type_fp ls1{p2,p3};
                  BOOST_TEST_INFO("intersection: " << bg::wkt(ls0) << " " << bg::wkt(ls1));
                  BOOST_CHECK_EQUAL(
                      path_finding::is_intersecting(p0, p1, p2, p3),
                      bg::intersects(ls0, ls1));
                }
              }
            }
          }
        }
      }
    }
  }
}

Answer 3

I have tried code it without testing... I hope it works! ^^

public static Point intersection(Segment s1, Segment s2) {
    // Components of the first segment's rect.
    Point v1 = new Point(s1.p2.x - s1.p1.x, s1.p2.y - s1.p1.y); // Directional vector
    double a1 = v.y;
    double b1 = -v.x;
    double c1 = v1.x * s1.p1.y - s1.v.y * s1.p1.x;

    // Components of the second segment's rect.
    Point v2 = new Point(s2.p2.x - s2.p1.x, s2.p2.y - s2.p1.y);
    double a2 = v2.y;
    double b2 = -v2.x;
    double c2 = v2.x * s2.p1.y - s2.v.y * s2.p1.x;

    // Calc intersection between RECTS.
    Point intersection = null;
    double det = a1 * b2 - b1 * a2;
    if (det != 0) {
        intersection = new Point(
            (b2 * (-c1) - b1 * (-c2)) / det;
            (a1 * (-c2) - a2 * (-c1)) / det;
        );
    }

    // Checks ff segments collides.
    if (
        intersection != null &&
        (
            (s1.p1.x <= intersection.x && intersection.x <= s1.p2.x) ||
            (s1.p2.x <= intersection.x && intersection.x <= s1.p1.x)
        ) &&
        (
            (s1.p1.y <= intersection.y && intersection.y <= s1.p2.y) ||
            (s1.p2.y <= intersection.y && intersection.y <= s1.p1.y)
        ) &&
        (
            (s2.p1.x <= intersection.x && intersection.x <= s2.p2.x) ||
            (s2.p2.x <= intersection.x && intersection.x <= s2.p1.x)
        ) &&
        (
            (s2.p1.y <= intersection.y && intersection.y <= s2.p2.y) ||
            (s2.p2.y <= intersection.y && intersection.y <= s2.p1.y)
        )
    )
        return intersection;

    return null;
};

Answer 4

Line intersection without special cases

Coming from a background of projective geometry, I'd write the points in homogeneous coordinates:

v1 = [x1, y1, 1]
v2 = [x2, y2, 1]
v3 = [x3, y3, 1]
v4 = [x4, y4, 1]

Then both the line joining two points and the intersection of two lines can be expressed using the cross product:

[x5, y5, z5] = (v1 × v2) × (v3 × v4)

which you can dehomogenize to find the resulting point as

[x5/z5, y5/z5]

without having to deal with any special cases. If your lines are parallel, the last point would lead to a division by zero, though, so you might want to catch that case.

Restriction to segments

The above is for infinite lines, though. You might want to keep the code which returns null if the point of intersection falls outside the bounding box. But if you want real segments, that code is incorrect: you could have a point of intersection which lies outside one of the segments but still inside the bounding box.

A proper check can be implemented using an orientation-checking predicate. The determinant of three of the vectors vi given above will have positive sign if the triangle they form has one orientation, and negative sign for the opposite orientation. So the points v3 and v4 lie on different sides of s1 if

det(v1, v2, v3) * det(v1, v2, v4) < 0

and in a similar way v1 and v2 lie on different sides of s2 if

det(v3, v4, v1) * det(v3, v4, v2) < 0

so if both of these are satisfied, you have an intersection between the segments. If you want to include the segment endpoints, change the < to a ≤ in these inequalities.