Как найти точку пересечения прямых python

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.

answered Mar 10, 2017 at 20:56

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

user1248490user1248490

9539 silver badges16 bronze badges

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

answered Aug 25, 2021 at 14:32

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!

answered Mar 31, 2017 at 17:36

dashesydashesy

2,5593 gold badges44 silver badges61 bronze badges

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:

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)

answered Sep 6, 2019 at 11:32

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)

answered Jan 22, 2022 at 17:13

TatarizeTatarize

10.1k4 gold badges57 silver badges64 bronze badges

xmduhanxmduhan

94712 silver badges14 bronze badges