# Numpy and line intersections

Posted on

### Question :

Numpy and line intersections

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.

``````#
# line segment intersection using vectors
# see Computer Graphics by F.S. Hill
#
from numpy import *
def perp( a ) :
b = empty_like(a)
b = -a
b = a
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, h)           # get first line
l2 = np.cross(h, h)           # 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 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.get_paths().vertices
B = Bcs.collections.get_paths().vertices

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

have fun

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
``````

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!

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 - p0
dx = q0 - p0
lhs0 = [-dy, dx]
rhs0 = p0 * dx - dy * p0
else:
lhs0 = [-m0, 1]
rhs0 = p0 - m0 * p0

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

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
``````

``````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)
``````