Python: ближайшая точка к линии

У меня такой вопрос. У меня есть поле с координатами и тремя точками, которые составляют линию. Теперь я хочу вычислить кратчайшее расстояние от всех координат коробки до этой линии. У меня есть три метода для этого, и версии vtk и numpy всегда дают один и тот же результат, но не метод расстояния shapely. Но мне нужна стройная версия, потому что я хочу измерить ближайшее расстояние от точки до всей линии, а не до отдельных сегментов линии. Вот пример кода. Итак, проблема в pdist:

from shapely.geometry import LineString, Point
import vtk, numpy as np
import itertools

import math

from numpy.linalg import norm

x1=np.arange(4,21)
y1=np.arange(4,21)
z1=np.arange(-7,6)

linepoints = np.array([[1,10,0],[10,10,0],[15,15,0]])


for i in itertools.product(x1,y1,z1):
    
    for m in range(len(linepoints)-1):
        
        line3 = LineString([linepoints[m],linepoints[m+1]])
        
        p = Point(i)
        
        d = norm(np.cross(linepoints[m]-linepoints[m+1], linepoints[m]-i))/norm(linepoints[m+1]-linepoints[m])
        
        dist=math.sqrt(vtk.vtkLine().DistanceToLine(i,linepoints[m],linepoints[m+1]))
        
        pdist = p.distance(line3)
        
        print(d,dist,pdist)

python numpy scipy shapely
person Varlor    schedule 23.05.2017    source источник

Ответы (1)


arrow_upward
3
arrow_downward

Проблема в том, что с помощью перекрестного произведения вы вычисляете ортогональное расстояние до линии, охватываемой сегментом, определенным точками linepoints[m] и linepoints[m+1]. С другой стороны, Shapely вычисляет расстояние до сегмента, то есть возвращает расстояние либо до ортогональной проекции, либо до одной из граничных точек, если ортогональная проекция выпадает «за пределы» сегмента.

Чтобы получить согласованные результаты, вы можете рассчитать ортогональную проекцию самостоятельно, а затем вызвать метод Shapely distance:

import numpy as np
from shapely.geometry import Point, LineString


A = np.array([1,0])
B = np.array([3,0])
C = np.array([0,1])


l = LineString([A, B])
p = Point(C)


d = np.linalg.norm(np.cross(B - A, C - A))/np.linalg.norm(B - A)

n = B - A
v = C - A

z = A + n*(np.dot(v, n)/np.dot(n, n))

print(l.distance(p), d, Point(z).distance(p))
#1.4142135623730951 1.0 1.0

Однако обратите внимание, что Shapely фактически игнорирует координату z. Так например:

import numpy as np
from shapely.geometry import Point, LineString

A = np.array([1,0,1])
B = np.array([0,0,0])

print(Point([1,0,1]).distance(Point([0,0,0])))

вернуть как расстояние просто 1.

РЕДАКТИРОВАТЬ: на основе вашего комментария здесь будет версия, которая вычисляет расстояние (для произвольного количества измерений) до сегмента:

from shapely.geometry import LineString, Point
import numpy as np
import itertools

import math

from numpy.linalg import norm

x1=np.arange(4,21)
y1=np.arange(4,21)
z1=np.arange(-7,6)

linepoints = np.array([[1,10,0],[10,10,0],[15,15,0]])

def dist(A, B, C):
    """Calculate the distance of point C to line segment spanned by points A, B.

    """

    a = np.asarray(A)
    b = np.asarray(B)
    c = np.asarray(C)

    #project c onto line spanned by a,b but consider the end points
    #should the projection fall "outside" of the segment    
    n, v = b - a, c - a

    #the projection q of c onto the infinite line defined by points a,b
    #can be parametrized as q = a + t*(b - a). In terms of dot-products,
    #the coefficient t is (c - a).(b - a)/( (b-a).(b-a) ). If we want
    #to restrict the "projected" point to belong to the finite segment
    #connecting points a and b, it's sufficient to "clip" it into
    #interval [0,1] - 0 corresponds to a, 1 corresponds to b.

    t = max(0, min(np.dot(v, n)/np.dot(n, n), 1))
    return np.linalg.norm(c - (a + t*n)) #or np.linalg.norm(v - t*n)


for coords in itertools.product(x1,y1,z1):
    for m in range(len(linepoints)-1):

        line3 = LineString([linepoints[m],linepoints[m+1]])
        d = dist(linepoints[m], linepoints[m+1], coords)
        print(coords, d)
person ewcz    schedule 23.05.2017
comment
Ах, это проблема! Что было бы тогда наилучшим подходом для расчета ближайшего расстояния от всех 3D-точек до полилинии? Потому что в некоторых случаях точки имеют самое близкое расстояние к бесконечному продолжению линейного сегмента, но мне нужно самое близкое расстояние до точки на сегменте / линии и, конечно же, в направлении z. Я думаю, это то, что вы описали с помощью той ортогональной проекции? - person Varlor; 23.05.2017
comment
@Varlor Я обновил ответ, используя, как я думаю, общий метод. В функции dist ограничение параметра t интервалом [0,1] гарантирует, что он вычисляет расстояние до конечного сегмента. Без этого ограничения он вернул бы расстояние до бесконечной линии ... - person ewcz; 23.05.2017
comment
Ах, здорово, работает. Можете ли вы объяснить, что именно происходит математически в строках: n, v = b - a, c - at = max (0, min (np.dot (v, n) /np.dot (n, n), 1)) return np.linalg.norm (c - (a + t n)) # или np.linalg.norm (v - t n) - person Varlor; 23.05.2017
comment
@Varlor Я добавил дополнительные комментарии, короче говоря, t определяет долю отрезка линии, соединяющего a и b, где расположена ортогональная проекция c. - person ewcz; 23.05.2017