Пересечения Numpy и Line

Как бы я использовал numpy для вычисления пересечения между двумя сегментами линии?

В коде у меня есть segment1 = ((x1,y1),(x2,y2)) и segment2 = ((x1,y1),(x2,y2)). Примечание segment1 не равно segment2. Итак, в моем коде я также вычислял наклон и y-перехват, было бы неплохо, если бы этого можно было избежать, но я не знаю, как это сделать.

Я использовал правило Крамера с функцией, которую я написал на Python, но я хотел бы найти более быстрый способ сделать это.


person Xavier    schedule 15.07.2010    source источник


Ответы (7)


Украдено непосредственно с сайта https://web.archive.org/web/20111108065352/https://www.cs.mun.ca/%7Erod/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)
person Hamish Grubijan    schedule 15.07.2010
comment
Спасибо за подсказку. Увидев этот алгоритм, я начал читать об этом. Вот хорошее объяснение IMO softsurfer.com/Archive/algorithm_0104/algorithm_0104B.htm . Надеюсь, это также послужит чьему-то любопытству. - person Maik Beckmann; 17.05.2011
comment
Примечание для тех, кто использует приведенный выше код: убедитесь, что вы передаете массив чисел с плавающей запятой в seg_intersect, иначе из-за округления могут произойти неожиданные вещи. - person schickm; 14.10.2012
comment
Кроме того, не забудьте проверить, равно ли denom нулю, иначе вы получите ошибку деления на ноль. (Это происходит, если прямые параллельны.) - person Gareth Rees; 21.12.2012
comment
@schickm, где возникает эта проблема с округлением? во время разделения? - person bennlich; 16.04.2015
comment
Справляется ли это с коллинеарностью? см.: stackoverflow. com/questions/3838329/ - person Nathan; 04.07.2018
comment
Предоставленная вами ссылка мертва (понятно, девять лет спустя...), но, к счастью, она была сохранена интернет-архивом. Это кажется полезным даже сейчас, поэтому вот ссылка на заархивированную версию: https://web.archive.org/web/20111108065352/https://www.cs.mun.ca/~rod/2500/notes/numpy-arrays/numpy-arrays.html - person cji; 13.03.2020

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

Объяснение

Обратите внимание, что уравнение линии ax+by+c=0. Итак, если точка находится на этой прямой, то это решение (a,b,c).(x,y,1)=0 (. — скалярное произведение)

пусть l1=(a1,b1,c1), l2=(a2,b2,c2) будут двумя линиями, а p1=(x1,y1,1), p2=(x2,y2,1) будут двумя точками.


Нахождение прямой, проходящей через две точки:

пусть t=p1xp2 (перекрестное произведение двух точек) будет вектором, представляющим линию.

Мы знаем, что p1 находится на линии t, потому что t.p1 = (p1xp2).p1=0. Мы также знаем, что p2 находится на t, потому что t.p2 = (p1xp2).p2=0. Таким образом, t должна быть линией, проходящей через p1 и p2.

Это означает, что мы можем получить векторное представление линии, взяв векторное произведение двух точек на этой линии.


Находим точку пересечения:

Теперь пусть r=l1xl2 (перекрестное произведение двух прямых) будет вектором, представляющим точку

Мы знаем, что r лежит на l1, потому что r.l1=(l1xl2).l1=0. Мы также знаем, что r лежит на l2, потому что r.l2=(l1xl2).l2=0. Таким образом, r должно быть точкой пересечения линий l1 и l2.

Интересно, что мы можем найти точку пересечения, взяв векторное произведение двух прямых.

person Norbu Tsering    schedule 10.03.2017
comment
Можете ли вы привести пример использования, начиная с 2 строк, каждая из которых определяется двумя 2D-точками? - person Matthias; 08.11.2017
comment
@Matthias Я добавил пример - person Norbu Tsering; 08.11.2017
comment
Спасибо! Но get_slope_intercept выдает исключение: одна линия горизонтальная, а другая вертикальная перпендикулярная. пример: (1, 1), (3, 1), (2,5, 2), (2,5, 0) - person Matthias; 08.11.2017
comment
Ах, это правильно. Вертикальные линии сделают матрицу коэффициентов единственной. Дай мне день. Я позабочусь об этом, когда у меня будет шанс - person Norbu Tsering; 08.11.2017
comment
@Matthias Обновил ответ. я думаю ты будешь очень доволен лол - person Norbu Tsering; 09.11.2017
comment
Какое красивое объяснение! - person Blanka; 06.03.2018
comment
Почему вы говорите, что t — это линия, проходящая через p1 и p2? Рассматривая эти точки как векторные смещения относительно начала координат в пространстве (начало координат равно [0,0], поэтому точка [x, y] является смещением от начала координат), когда вы берете векторное произведение между этими двумя векторами смещения вы получаете другой вектор, параллельный им и выходящий за пределы экрана, который не является вектором, проходящим через точки. - person R. Navega; 08.01.2019
comment
@Р. Навега, потому что t.p1=t.p2=0 - person Norbu Tsering; 08.08.2020
comment
На самом деле мы работаем только в одной плоскости трехмерного пространства. Плоскость состоит из всех точек с третьей составляющей, равной 1. Вот почему я говорю, что это вектор, представляющий линию. - person Norbu Tsering; 08.08.2020
comment
@NorbuTsering хорошо, поправьте меня, если я ошибаюсь, t на самом деле является вектором коэффициентов линии, таким же, как l1 и l2, которые вы используете. Эти векторы не могут быть визуализированы как линии сами по себе, вам нужно визуализировать их, используя точки [x, y], которые решают их линейное уравнение, например, для t это будут точки, которые решают t₁ * x + t₂ * y + t₃ = 0 (где t является 3- компонентный вектор) - person R. Navega; 09.08.2020
comment
Правильный. Iirc, я получил это из книги о проективной геометрии. - person Norbu Tsering; 10.08.2020

Возможно, это поздний ответ, но он был первым хитом, когда я погуглил «пересечения numpy line». В моем случае у меня есть две линии на плоскости, и я хотел быстро получить любые пересечения между ними, а решение Хэмиша было бы медленным — требовалось бы вложенный цикл for для всех сегментов линии.

Вот как это сделать без цикла for (это довольно быстро):

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

Затем, чтобы использовать его, укажите две строки в качестве аргументов, где arg — это матрица с двумя столбцами, каждая строка соответствует точке (x, y):

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

развлекайся

person marmaduke    schedule 02.02.2012
comment
переменная m1inv не используется, это нормально? - person adrienlucca.net; 29.12.2014
comment
Что вы подразумеваете под любыми пересечениями между ними? сколько их там? - person Gulzar; 27.12.2018

Это версия ответа @Hamish Grubijan, которая также работает для нескольких точек в каждом из входных аргументов, т. е. a1, a2, b1, b2 могут быть массивами строк Nx2 из 2D-точек. Функция perp заменяется скалярным произведением.

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
person user1248490    schedule 16.11.2016

Вот (немного вынужденный) однострочник:

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)

Интерполируйте разницу (по умолчанию линейная) и найдите 0 обратной.

Ваше здоровье!

person Andy Reagan    schedule 31.03.2017
comment
Извините за комментарий к старому сообщению, но как это должно работать? Вы не можете вычесть кортежи, и использование массивов np возвращает ошибку, что x (segment1) не может иметь несколько измерений. - person justry; 23.04.2020
comment
Да хороший вопрос. Мне пришлось действительно подумать, я отредактировал свой ответ, включив данные. Короче говоря, это просто возвращает значение x. - person Andy Reagan; 24.04.2020
comment
Я не уверен, как это будет работать для двух сегментов с отдельными координатами x и Y, но у меня это сработало, поскольку все, что я хотел, это пересечение с базовой линией. Спасибо! - person justry; 25.04.2020

Это то, что я использую для поиска пересечения линий, оно работает либо с двумя точками каждой линии, либо просто с точкой и ее наклоном. Я в основном решаю систему линейных уравнений.

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
person dashesy    schedule 17.10.2014

person    schedule
comment
Перекрестное произведение Numpy слишком медленное. это занимает 47 секунд, а мое решение - 800 мс. - person Sadekujjaman; 11.09.2019
comment
Это было бы более полезно с каким-то объяснением того, чем это отличается от других ответов. - person RTbecard; 19.01.2020