Numpy and line intersections

Question

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 segment 1 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.

Answer 1

Stolen directly from http://www.cs.mun.ca/~rod/2500/notes/numpy-arrays/numpy-arrays.html

#
# line segment intersection using vectors
# see Computer Graphics by F.S. Hill
#
from numpy import *
def perp( a ) :
    b = empty_like(a)
    b[0] = -a[1]
    b[1] = a[0]
    return b

# line segment a given by endpoints a1, a2
# line segment b given by endpoints b1, b2
# return 
def seg_intersect(a1,a2, b1,b2) :
    da = a2-a1
    db = b2-b1
    dp = a1-b1
    dap = perp(da)
    denom = dot( dap, db)
    num = dot( dap, dp )
    return (num / denom.astype(float))*db + b1

p1 = array( [0.0, 0.0] )
p2 = array( [1.0, 0.0] )

p3 = array( [4.0, -5.0] )
p4 = array( [4.0, 2.0] )

print seg_intersect( p1,p2, p3,p4)

p1 = array( [2.0, 2.0] )
p2 = array( [4.0, 3.0] )

p3 = array( [6.0, 0.0] )
p4 = array( [6.0, 3.0] )

print seg_intersect( p1,p2, p3,p4)

Answer 2

import numpy as np

def get_intersect(a1, a2, b1, b2):
    """ 
    Returns the point of intersection of the lines passing through a2,a1 and b2,b1.
    a1: [x, y] a point on the first line
    a2: [x, y] another point on the first line
    b1: [x, y] a point on the second line
    b2: [x, y] another point on the second line
    """
    s = np.vstack([a1,a2,b1,b2])        # s for stacked
    h = np.hstack((s, np.ones((4, 1)))) # h for homogeneous
    l1 = np.cross(h[0], h[1])           # get first line
    l2 = np.cross(h[2], h[3])           # get second line
    x, y, z = np.cross(l1, l2)          # point of intersection
    if z == 0:                          # lines are parallel
        return (float('inf'), float('inf'))
    return (x/z, y/z)

if __name__ == "__main__":
    print get_intersect((0, 1), (0, 2), (1, 10), (1, 9))  # parallel  lines
    print get_intersect((0, 1), (0, 2), (1, 10), (2, 10)) # vertical and horizontal lines
    print get_intersect((0, 1), (1, 2), (0, 10), (1, 9))  # another line for fun

---

Explanation

Note that the equation of a line is ax+by+c=0. So if a point is on this line, then it is a solution to (a,b,c).(x,y,1)=0 (. is the dot product)

let l1=(a1,b1,c1), l2=(a2,b2,c2) be two lines and p1=(x1,y1,1), p2=(x2,y2,1) be two points.

Finding the line passing through two points:

let t=p1xp2 (the cross product of two points) be a vector representing a line.

We know that p1 is on the line t because t.p1 = (p1xp2).p1=0. We also know that p2 is on t because t.p2 = (p1xp2).p2=0. So t must be the line passing through p1 and p2.

This means that we can get the vector representation of a line by taking the cross product of two points on that line.

Finding the point of intersection:

Now let r=l1xl2 (the cross product of two lines) be a vector representing a point

We know r lies on l1 because r.l1=(l1xl2).l1=0. We also know r lies on l2 because r.l2=(l1xl2).l2=0. So r must be the point of intersection of the lines l1 and l2.

Interestingly, we can find the point of intersection by taking the cross product of two lines.

Answer 3

This is is a late response, perhaps, but it was the first hit when I Googled 'numpy line intersections'. In my case, I have two lines in a plane, and I wanted to quickly get any intersections between them, and Hamish's solution would be slow -- requiring a nested for loop over all line segments.

Here's how to do it without a for loop (it's quite fast):

from numpy import where, dstack, diff, meshgrid

def find_intersections(A, B):

    # min, max and all for arrays
    amin = lambda x1, x2: where(x1<x2, x1, x2)
    amax = lambda x1, x2: where(x1>x2, x1, x2)
    aall = lambda abools: dstack(abools).all(axis=2)
    slope = lambda line: (lambda d: d[:,1]/d[:,0])(diff(line, axis=0))

    x11, x21 = meshgrid(A[:-1, 0], B[:-1, 0])
    x12, x22 = meshgrid(A[1:, 0], B[1:, 0])
    y11, y21 = meshgrid(A[:-1, 1], B[:-1, 1])
    y12, y22 = meshgrid(A[1:, 1], B[1:, 1])

    m1, m2 = meshgrid(slope(A), slope(B))
    m1inv, m2inv = 1/m1, 1/m2

    yi = (m1*(x21-x11-m2inv*y21) + y11)/(1 - m1*m2inv)
    xi = (yi - y21)*m2inv + x21

    xconds = (amin(x11, x12) < xi, xi <= amax(x11, x12), 
              amin(x21, x22) < xi, xi <= amax(x21, x22) )
    yconds = (amin(y11, y12) < yi, yi <= amax(y11, y12),
              amin(y21, y22) < yi, yi <= amax(y21, y22) )

    return xi[aall(xconds)], yi[aall(yconds)]

Then to use it, provide two lines as arguments, where is arg is a 2 column matrix, each row corresponding to an (x, y) point:

# example from matplotlib contour plots
Acs = contour(...)
Bsc = contour(...)

# A and B are the two lines, each is a 
# two column matrix
A = Acs.collections[0].get_paths()[0].vertices
B = Bcs.collections[0].get_paths()[0].vertices

# do it
x, y = find_intersections(A, B)

have fun

Answer 4

This is a version of @Hamish Grubijan's answer that also works for multiple points in each of the input arguments, i.e., a1, a2, b1, b2 can be Nx2 row arrays of 2D points. The perp function is replaced by a dot product.

T = np.array([[0, -1], [1, 0]])
def line_intersect(a1, a2, b1, b2):
    da = np.atleast_2d(a2 - a1)
    db = np.atleast_2d(b2 - b1)
    dp = np.atleast_2d(a1 - b1)
    dap = np.dot(da, T)
    denom = np.sum(dap * db, axis=1)
    num = np.sum(dap * dp, axis=1)
    return np.atleast_2d(num / denom).T * db + b1

Answer 5

Here's a (bit forced) one-liner:

import numpy as np
from scipy.interpolate import interp1d

x = np.array([0, 1])
segment1 = np.array([0, 1])
segment2 = np.array([-1, 2])

x_intersection = interp1d(segment1 - segment2, x)(0)
# if you need it:
y_intersection = interp1d(x, segment1)(x_intersection)

Interpolate the difference (default is linear), and find a 0 of the inverse.

Cheers!

Answer 6

This is what I use to find line intersection, it works having either 2 points of each line, or just a point and its slope. I basically solve the system of linear equations.

def line_intersect(p0, p1, m0=None, m1=None, q0=None, q1=None):
    ''' intersect 2 lines given 2 points and (either associated slopes or one extra point)
    Inputs:
        p0 - first point of first line [x,y]
        p1 - fist point of second line [x,y]
        m0 - slope of first line
        m1 - slope of second line
        q0 - second point of first line [x,y]
        q1 - second point of second line [x,y]
    '''
    if m0 is  None:
        if q0 is None:
            raise ValueError('either m0 or q0 is needed')
        dy = q0[1] - p0[1]
        dx = q0[0] - p0[0]
        lhs0 = [-dy, dx]
        rhs0 = p0[1] * dx - dy * p0[0]
    else:
        lhs0 = [-m0, 1]
        rhs0 = p0[1] - m0 * p0[0]

    if m1 is  None:
        if q1 is None:
            raise ValueError('either m1 or q1 is needed')
        dy = q1[1] - p1[1]
        dx = q1[0] - p1[0]
        lhs1 = [-dy, dx]
        rhs1 = p1[1] * dx - dy * p1[0]
    else:
        lhs1 = [-m1, 1]
        rhs1 = p1[1] - m1 * p1[0]

    a = np.array([lhs0, 
                  lhs1])

    b = np.array([rhs0, 
                  rhs1])
    try:
        px = np.linalg.solve(a, b)
    except:
        px = np.array([np.nan, np.nan])

    return px

Answer 7

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)