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
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
. is the dot product)
l2=(a2,b2,c2) be two lines and
p2=(x2,y2,1) be two points.
Finding the line passing through two points:
t=p1xp2 (the cross product of two points) be a vector representing a line.
We know that
p1 is on the line
t.p1 = (p1xp2).p1=0.
We also know that
p2 is on
t.p2 = (p1xp2).p2=0. So
t must be the line passing through
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:
r=l1xl2 (the cross product of two lines) be a vector representing a point
r lies on
r.l1=(l1xl2).l1=0. We also know
r lies on
r must be the point of intersection of the lines
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)
This is a version of @Hamish Grubijan’s answer that also works for multiple points in each of the input arguments, i.e.,
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.
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)