Numpy and line intersections

Viewed 83058

How would I use numpy to calculate the intersection between two line segments?

In the code I have segment1 = ((x1,y1),(x2,y2)) and segment2 = ((x1,y1),(x2,y2)). Note segment1 does not equal segment2. So in my code I've also been calculating the slope and y-intercept, it would be nice if that could be avoided but I don't know of a way how.

I've been using Cramer's rule with a function I wrote up in Python but I'd like to find a faster way of doing this.

10 Answers

I would like to add something small here. The original question is about line segments. I arrived here, because I was looking for line segment intersection, which in my case meant that I need to filter those cases, where no intersection of the line segments exists. Here is some code which does that:

def line_intersection(x1, y1, x2, y2, x3, y3, x4, y4):
    """find the intersection of line segments A=(x1,y1)/(x2,y2) and
    B=(x3,y3)/(x4,y4). Returns a point or None"""
    denom = ((x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4))
    if denom==0: return None
    px = ((x1 * y2 - y1 * x2) * (x3 - x4) - (x1 - x2) * (x3 * y4 - y3 * x4)) / denom
    py = ((x1 * y2 - y1 * x2) * (y3 - y4) - (y1 - y2) * (x3 * y4 - y3 * x4)) / denom
    if (px - x1) * (px - x2) < 0 and (py - y1) * (py - y2) < 0 \
      and (px - x3) * (px - x4) < 0 and (py - y3) * (py - y4) < 0:
        return [px, py]
    else:
        return None

In case you are looking for a vectorized version where we can rule out vertical line segments.

def intersect(a):
    # a numpy array with dimension [n, 2, 2, 2]
    # axis 0: line-pair, axis 1: two lines, axis 2: line delimiters axis 3: x and y coords
    # for each of the n line pairs a boolean is returned stating of the two lines intersect
    # Note: the edge case of a vertical line is not handled.
    m = (a[:, :, 1, 1] - a[:, :, 0, 1]) / (a[:, :, 1, 0] - a[:, :, 0, 0])
    t = a[:, :, 0, 1] - m[:, :] * a[:, :, 0, 0]
    x = (t[:, 0] - t[:, 1]) / (m[:, 1] - m[:, 0])
    y = m[:, 0] * x + t[:, 0]
    r = a.min(axis=2).max(axis=1), a.max(axis=2).min(axis=1)
    return (x >= r[0][:, 0]) & (x <= r[1][:, 0]) & (y >= r[0][:, 1]) & (y <= r[1][:, 1])

A sample invocation would be:

intersect(np.array([
    [[[1, 2], [2, 2]],
     [[1, 2], [1, 1]]], # I
    [[[3, 4], [4, 4]],
     [[4, 4], [5, 6]]], # II
    [[[2, 0], [3, 1]],
     [[3, 0], [4, 1]]], # III
    [[[0, 5], [2, 5]],
     [[2, 4], [1, 3]]], # IV
]))
# returns [False, True, False, False]

Visualization (I need more reputation to post images here).

We can solve this 2D line intersection problem using determinant.

To solve this, we have to convert our lines to the following form: ax+by=c. where

 a = y1 - y2
 b = x1 - x2
 c = ax1 + by1 

If we apply this equation for each line, we will got two line equation. a1x+b1y=c1 and a2x+b2y=c2.

Now when we got the expression for both lines.
First of all we have to check if the lines are parallel or not. To examine this we want to find the determinant. The lines are parallel if the determinant is equal to zero.
We find the determinant by solving the following expression:

det = a1 * b2 - a2 * b1

If the determinant is equal to zero, then the lines are parallel and will never intersect. If the lines are not parallel, they must intersect at some point.
The point of the lines intersects are found using the following formula:

equation for finding line intersection using determinant

class Point:
    def __init__(self, x, y):
        self.x = x
        self.y = y


'''
finding intersect point of line AB and CD 
where A is the first point of line AB
and B is the second point of line AB
and C is the first point of line CD
and D is the second point of line CD
'''



def get_intersect(A, B, C, D):
    # a1x + b1y = c1
    a1 = B.y - A.y
    b1 = A.x - B.x
    c1 = a1 * (A.x) + b1 * (A.y)

    # a2x + b2y = c2
    a2 = D.y - C.y
    b2 = C.x - D.x
    c2 = a2 * (C.x) + b2 * (C.y)

    # determinant
    det = a1 * b2 - a2 * b1

    # parallel line
    if det == 0:
        return (float('inf'), float('inf'))

    # intersect point(x,y)
    x = ((b2 * c1) - (b1 * c2)) / det
    y = ((a1 * c2) - (a2 * c1)) / det
    return (x, y)

I wrote a module for line to compute this and some other simple line operations. It is implemented in c++, so it works very fast. You can install FastLine via pip and then use it in this way:

from FastLine import Line
# define a line by two points
l1 = Line(p1=(0,0), p2=(10,10))
# or define a line by slope and intercept
l2 = Line(m=0.5, b=-1)

# compute intersection
p = l1.intersection(l2)
# returns (-2.0, -2.0)
Related