Существует ли быстрый Python алгоритм для нахождения всех точек в наборе данных, которые l ie в данном круге?

Question

У меня большой объем данных (более 10 ^ 5 баллов). Я ищу быстрый алгоритм, который находит все точки в наборе данных, которые l ie в круге, заданном его центральной точкой и радиусом.

Я думал об использовании дерева kd для вычисления, например, 10 Ближайшие точки к центру круга, а затем проверьте, находятся ли они внутри круга. Но я не уверен, что это правильный путь.

Answer 1

Я сравнил версию numexpr с простой реализацией Numpy следующим образом:

#!/usr/bin/env python3

import numpy as np
import numexpr as ne

# Ensure repeatable, deterministic randomness!
np.random.seed(42)

# Generate test arrays
N = 1000000
X = np.random.rand(N)
Y = np.random.rand(N)

# Define centre and radius
cx = cy = r = 0.5

def method1(X,Y,cx,cy,r):
    """Straight Numpy determination of points in circle"""
    d = (X-cx)**2 + (Y-cy)**2
    res = d < r**2
    return res

def method2(X,Y,cx,cy,r):
    """Numexpr determination of points in circle"""
    res = ne.evaluate('((X-cx)**2 + (Y-cy)**2)<r**2')
    return res

def method3(data,a,b,r): 
    """List based determination of points in circle, with pre-filtering using a square"""
    in_square_points = [(x,y) for (x,y) in data if a-r < x < a+r and b-r < y < b+r] 
    in_circle_points = [(x,y) for (x,y) in in_square_points if (x-a)**2 + (y-b)**2 < r**2] 
    return in_circle_points

# Timing
%timeit method1(X,Y,cx,cy,r)

%timeit method2(X,Y,cx,cy,r)

# Massage input data (before timing) to match agorithm
data=[(x,y) for x,y in zip(X,Y)]
%timeit method3(data,cx,cy,r)

Затем я рассчитал ее в I Python следующим образом:

%timeit method1(X,Y,cx,cy,r)                                                                                                                     
6.68 ms ± 246 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

%timeit method2(X,Y,cx,cy,r)                                                                                                                     
743 µs ± 17.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

%timeit method3(data,cx,cy,r)
1.11 s ± 9.81 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Таким образом, версия numexpr вышла в 9 раз быстрее. Поскольку точки l ie находятся в диапазоне [0..1], алгоритм эффективно вычисляет pi, и оба метода получаются одинаковыми:

method1(X,Y,cx,cy,r).sum()                                                                                                                       
784973

method2(X,Y,cx,cy,r).sum()                                                                                                                       
784973

len(method3(data,cx,cy,r))                                                                                                                       
784973

4 * 784973 / N
3.139

Примечание : Я должен отметить, что numexpr автоматически выполняет многопоточность вашего кода на нескольких ядрах процессора. Если вам хочется поэкспериментировать с количеством потоков, вы можете изменить его динамически перед вызовом method2() или даже внутри него:

# Split calculations across 6 threads
ne.set_num_threads(6) 

---

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

Answer 2

Если вас интересует только количество точек в кружке, вы можете попробовать Numba.

import numpy as np
import numba as nb
import numexpr as ne

def method2(X,Y,cx,cy,r):
    """Numexpr method"""
    res = ne.evaluate('((X-cx)**2 + (Y-cy)**2) < r**2')
    return res

@nb.njit(fastmath=True,parallel=True)
def method3(X,Y,cx,cy,r):
    acc=0
    for i in nb.prange(X.shape[0]):
        if ((X[i]-cx)**2 + (Y[i]-cy)**2) < r**2:
            acc+=1
    return acc

Время

# Ensure repeatable, deterministic randomness!
np.random.seed(42)

# Generate test arrays
N = 1000000
X = np.random.rand(N)
Y = np.random.rand(N)

# Define centre and radius
cx = cy = r = 0.5

#@Mark Setchell
%timeit method2(X,Y,cx,cy,r)
#825 µs ± 22.1 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit method3(X,Y,cx,cy,r)
#480 µs ± 94.4 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

Answer 3

Структура данных с разделением пространства, такая как дерево kd или quadtree, может точно ответить на ваш запрос; нет необходимости в heuristi c, например, "10 ближайших точек".

Вы можете рекурсивно искать дерево:

1. Если ограничивающий прямоугольник текущего узла вообще не пересекает окружность , игнорируйте его.
2. В противном случае, если ограничивающий прямоугольник текущего узла полностью содержится внутри круга, вернуть все точки в нем.
3. В противном случае, если ограничивающий прямоугольник текущего узла частично перекрывается с окружностью :
   • Если текущий узел является листом, проверьте каждую точку по отдельности, чтобы увидеть, находится ли она внутри круга, и верните те, которые находятся. вернуть все точки, возвращаемые рекурсивными вызовами.

На шаге 2 прямоугольник полностью содержится в круге тогда и только тогда, когда его четыре угла равны. Для шагов 1 и 3 вы можете либо проверить, пересекает ли прямоугольник ограничивающую рамку круга, и консервативно выполнить рекурсию, потому что прямоугольник может пересекать круг; или вы можете сделать точную (но немного медленнее) проверку между прямоугольником и фактическим кругом.

Поскольку это Python, вместо того, чтобы возвращать точки в коллекции и затем объединять Объединяя результаты рекурсивных вызовов, вы можете получить очки из функции рекурсивного генератора. Это несколько упрощает реализацию.

Answer 4

Если вы хотите сначала отфильтровать большое количество вашего набора данных без значительных вычислений, вы можете использовать квадрат размера (2r x 2r) с тем же центром, что и круг (где r - радиус круга).

Посмотрите на эту картинку:

Если у вас есть координаты центра (a, b) и r радиуса, то точки (x, y) внутри квадрат подтвердите:

in_square_points = [(x,y) for (x,y) in data if a-r < x < a+r and b-r < y < b+r]

И, наконец, после этого фильтра вы можете применить уравнение окружности:

in_circle_points = [(x,y) for (x,y) in in_square_points if (x-a)**2 + (y-b)**2 < r**2]

** РЕДАКТИРОВАТЬ **

если ваш вход структурирован так:

data = [
    [13, 45],
    [-1, 2],
    ...
    [60, -4]
]

Тогда вы можете попробовать, если предпочитаете обычные циклы for:

in_square_points = []
for i in range(len(data)):
    x = data[i][0]
    y = data[i][1]
    if a-r < x < a+r and b-r < y < b+r:
        in_square_points.append([x, y])
print(in_square_points)

Answer 5

Чтобы проверить, находится ли точка (a, b) внутри круга с центром (x, y) и радиусом r, вы можете просто выполнить вычисление:

within_circle = ((x-a)**2 + (y-b)**2) <= r*r)

В этом уравнении используется свойство окружность, по которой можно получить абсолютное расстояние до точки (которая также используется в формуле расстояния, если вы заметили).