Здесь у вас есть решение, использующее rtree
библиотеку. Идея состоит в том, чтобы построить блоки, содержащие сегменты по диагонали, и использовать эти блоки для построения rtree. Это будет самая затратная по времени операция. Позже вы запрашиваете rtree с прямоугольником, центрированным в точке. Вы получаете несколько совпадений, которые вам нужно проверить по минимуму, но количество совпадений будет (надеюсь) на несколько порядков меньше, чем проверка по всем сегментам.
В solutions
dict вы получите для каждой точки идентификатор линии, ближайший сегмент, ближайшую точку (точку сегмента) и расстояние до точки.
В коде есть несколько комментариев, которые могут вам помочь. Учтите, что вы можете сериализовать rtree для дальнейшего использования. Фактически, я бы рекомендовал построить rtree, сохранить его, а затем использовать. Потому что исключения для корректировок констант MIN_SIZE
и INFTY
, вероятно, увеличатся, и вы не захотите потерять все вычисления, которые вы выполнили при построении rtree.
Слишком маленький MIN_SIZE
будет означать, что у вас могут быть ошибки в решениях, потому что, если прямоугольник вокруг точки не пересекает сегмент, он может пересекать прямоугольник сегмента, который не является ближайшим сегментом (легко придумать случай).
Слишком большой MIN_SIZE
будет означать слишком много ложных срабатываний, что в крайнем случае заставит код попробовать со всеми сегментами, и вы будете в том же положении, что и раньше, или в худшем случае, потому что теперь вы строите rtree вы действительно не используете.
Если данные являются реальными данными из города, я полагаю, вы знаете, что любой адрес будет пересекать сегмент с расстоянием меньше нескольких кварталов. Это сделает поиск практически логарифмическим.
Еще один комментарий. Я предполагаю, что нет слишком больших сегментов. Поскольку мы используем сегменты в качестве диагоналей прямоугольников в rtree, если у вас есть несколько больших сегментов в строке, это будет означать, что этому сегменту будет назначен огромный прямоугольник, и все поля адресов будут пересекать его. Чтобы избежать этого, вы всегда можете искусственно увеличить разрешение LineStrins, добавив больше промежуточных точек.
import math
from rtree import index
from shapely.geometry import Polygon, LineString
INFTY = 1000000
MIN_SIZE = .8
# MIN_SIZE should be a vaule such that if you build a box centered in each
# point with edges of size 2*MIN_SIZE, you know a priori that at least one
# segment is intersected with the box. Otherwise, you could get an inexact
# solution, there is an exception checking this, though.
def distance(a, b):
return math.sqrt( (a[0]-b[0])**2 + (a[1]-b[1])**2 )
def get_distance(apoint, segment):
a = apoint
b, c = segment
# t = <a-b, c-b>/|c-b|**2
# because p(a) = t*(c-b)+b is the ortogonal projection of vector a
# over the rectline that includes the points b and c.
t = (a[0]-b[0])*(c[0]-b[0]) + (a[1]-b[1])*(c[1]-b[1])
t = t / ( (c[0]-b[0])**2 + (c[1]-b[1])**2 )
# Only if t 0 <= t <= 1 the projection is in the interior of
# segment b-c, and it is the point that minimize the distance
# (by pitagoras theorem).
if 0 < t < 1:
pcoords = (t*(c[0]-b[0])+b[0], t*(c[1]-b[1])+b[1])
dmin = distance(a, pcoords)
return pcoords, dmin
elif t <= 0:
return b, distance(a, b)
elif 1 <= t:
return c, distance(a, c)
def get_rtree(lines):
def generate_items():
sindx = 0
for lid, l in lines:
for i in xrange(len(l)-1):
a, b = l[i]
c, d = l[i+1]
segment = ((a,b), (c,d))
box = (min(a, c), min(b,d), max(a, c), max(b,d))
#box = left, bottom, right, top
yield (sindx, box, (lid, segment))
sindx += 1
return index.Index(generate_items())
def get_solution(idx, points):
result = {}
for p in points:
pbox = (p[0]-MIN_SIZE, p[1]-MIN_SIZE, p[0]+MIN_SIZE, p[1]+MIN_SIZE)
hits = idx.intersection(pbox, objects='raw')
d = INFTY
s = None
for h in hits:
nearest_p, new_d = get_distance(p, h[1])
if d >= new_d:
d = new_d
s = (h[0], h[1], nearest_p, new_d)
result[p] = s
print s
#some checking you could remove after you adjust the constants
if s == None:
raise Exception("It seems INFTY is not big enough.")
pboxpol = ( (pbox[0], pbox[1]), (pbox[2], pbox[1]),
(pbox[2], pbox[3]), (pbox[0], pbox[3]) )
if not Polygon(pboxpol).intersects(LineString(s[1])):
msg = "It seems MIN_SIZE is not big enough. "
msg += "You could get inexact solutions if remove this exception."
raise Exception(msg)
return result
Я тестировал функции на этом примере.
xcoords = [i*10.0/float(1000) for i in xrange(1000)]
l1 = [(x, math.sin(x)) for x in xcoords]
l2 = [(x, math.cos(x)) for x in xcoords]
points = [(i*10.0/float(50), 0.8) for i in xrange(50)]
lines = [('l1', l1), ('l2', l2)]
idx = get_rtree(lines)
solutions = get_solution(idx, points)
И получил:
person
eguaio
schedule
20.09.2017