Numpy and line intersections

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

Stolen directly from https://web.archive.org/web/20111108065352/https://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)

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.

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