2016-12-24 5 views
2

Для списка точек (x, y), я пытаюсь найти ближайшие точки для каждой точки.Как индексировать список точек для более быстрого поиска близлежащих точек?

from collections import defaultdict 
from math import sqrt 
from random import randint 

# Generate a list of random (x, y) points 
points = [(randint(0, 100), randint(0, 100)) for _ in range(1000)] 

def is_nearby(point_a, point_b, max_distance=5): 
    """Two points are nearby if their Euclidean distance is less than max_distance""" 
    distance = sqrt((point_b[0] - point_a[0])**2 + (point_b[1] - point_a[1])**2) 
    return distance < max_distance 

# For each point, find nearby points that are within a radius of 5 
nearby_points = defaultdict(list) 
for point in points: 
    for neighbour in points: 
     if point != neighbour: 
      if is_nearby(point, neighbour): 
       nearby_points[point].append(neighbour) 

Есть ли способ, что я могу индекс points сделать выше поиск быстрее? Я считаю, что должен быть какой-то более быстрый способ, чем O (len (points) ** 2).

Edit: точки в целом может быть поплавки, а не только INTS

+0

Если ваша сетка всего 100 * 100, вы можете расположить свои точки внутри этой сетки. таким образом, вы можете резко сократить пространство поиска. –

+0

http://gis.stackexchange.com/questions/22082/how-can-i-use-r-tree-to-find-points-within-a-distance-in-spatialite –

ответ

1

это версия с фиксированной сеткой, где каждый Gridpoint удерживает количество образцов, которые находятся там.

поиск может быть сведен к простому пространству вокруг рассматриваемой точки.

from random import randint 
import math 

N = 100 
N_SAMPLES = 1000 

# create the grid 
grd = [[0 for _ in range(N)] for __ in range(N)] 

# set the number of points at a given gridpoint 
for _ in range(N_SAMPLES): 
    grd[randint(0, 99)][randint(0, 99)] += 1 

def find_neighbours(grid, point, distance): 

    # this will be: (x, y): number of points there 
    points = {} 

    for x in range(point[0]-distance, point[0]+distance): 
     if x < 0 or x > N-1: 
      continue 
     for y in range(point[1]-distance, point[1]+distance): 
      if y < 0 or y > N-1: 
       continue 
      dst = math.hypot(point[0]-x, point[1]-y) 
      if dst > distance: 
       continue 
      if grd[x][y] > 0: 
       points[(x, y)] = grd[x][y] 
    return points 

print(find_neighbours(grid=grd, point=(45, 36), distance=5)) 
# -> {(44, 37): 1, (45, 33): 1, ...} 
# meadning: there is one neighbour at (44, 37) etc... 

для дальнейшего optimzation: тесты для x и y могут быть предварительно вычислены для данного gridsize - math.hypot(point[0]-x, point[1]-y) не должно быть сделано, то для каждой точки.

, и это может быть хорошей идеей заменить сетку массивом numpy.


UPDATE

, если ваши точки float s вы можете создать int сетку, чтобы уменьшить пространство поиска:

from random import uniform 
from collections import defaultdict 
import math 

class Point: 
    def __init__(self, x, y): 
     self.x = x 
     self.y = y 

    @property 
    def x_int(self): 
     return int(self.x) 

    @property 
    def y_int(self): 
     return int(self.y) 

    def __str__(self): 
     fmt = '''{0.__class__.__name__}(x={0.x:5.2f}, y={0.y:5.2f})''' 
     return fmt.format(self) 

N = 100 
MIN = 0 
MAX = N-1 

N_SAMPLES = 1000 


# create the grid 
grd = [[[] for _ in range(N)] for __ in range(N)] 

# set the number of points at a given gridpoint 
for _ in range(N_SAMPLES): 
    p = Point(x=uniform(MIN, MAX), y=uniform(MIN, MAX)) 
    grd[p.x_int][p.y_int].append(p) 


def find_neighbours(grid, point, distance): 

    # this will be: (x_int, y_int): list of points 
    points = defaultdict(list) 

    # need to cast a slightly bigger net on the upper end of the range; 
    # int() rounds down 
    for x in range(point[0]-distance, point[0]+distance+1): 
     if x < 0 or x > N-1: 
      continue 
     for y in range(point[1]-distance, point[1]+distance+1): 
      if y < 0 or y > N-1: 
       continue 
      dst = math.hypot(point[0]-x, point[1]-y) 
      if dst > distance + 1: # account for rounding... is +1 enough? 
       continue 
      for pt in grd[x][y]: 
       if math.hypot(pt.x-x, pt.y-y) <= distance: 
        points[(x, y)].append(pt) 
    return points 

res = find_neighbours(grid=grd, point=(45, 36), distance=5) 

for int_point, points in res.items(): 
    print(int_point) 
    for point in points: 
     print(' ', point) 

выход выглядит примерно так:

(44, 36) 
    Point(x=44.03, y=36.93) 
(41, 36) 
    Point(x=41.91, y=36.55) 
    Point(x=41.73, y=36.53) 
    Point(x=41.56, y=36.88) 
... 

для удобства Points - теперь класс. может не быть необходимыми, хотя ...

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

также функцию find_neighbours принимает начальную point, состоящий из int s только в этой версии. это также может быть уточнено.

и есть много возможностей для улучшения: ось y может быть ограничена использованием тригонометрии. и для точек внутри круга нет необходимости в индивидуальной проверке; детальная проверка должна выполняться только близко к внешнему краю круга.

+0

Спасибо, а что, если точки являются float вместо ints? Этот метод работал бы только в том случае, если бы мы объединили float в ints – mchen

+0

Я думаю, что настройка вашего вышеуказанного подхода будет работать. Вместо поиска по фиксированной сетке поиск между 'bisect_left ((точка [0] +/- расстояние, точка [1] +/- расстояние), точки)' – mchen