Эффективный метод подсчета количества точек данных внутри сферы фиксированного радиуса с центром в каждой точке данных

Question

У меня есть база данных со многими точками данных, каждая из которых имеет координаты x, y, z.Я хочу посчитать количество точек, которые находятся на определенном расстоянии от соседних точек.Некоторые точки будут иметь пару в радиусе R, другие не будут.Я просто хочу посчитать количество пар на некотором расстоянии.Я мог бы легко написать алгоритм, чтобы сделать это, но он не был бы достаточно эффективным (поскольку я перебрал бы каждую точку данных).

Это похоже на то, что уже должно существовать в астропии, скипе и т. Д., Но я не могу найти то, что ищу.Есть ли что-нибудь, что достигает этого?

Answer 1

Как отметил @Davis Herring в комментариях, эффективный вариант - это дерево k-d.

k-d дерево - это алгоритм, который избегает подхода грубой силы и позволяет эффективно вычислять расстояния * (см. Нижнюю часть ответа для фона).

Существует несколько реализаций Python, одна из которых через SciPy:

Дерево SciPy k-d в Cython (быстрее, так как использует C / Cython)

Дерево SciPy k-d в чистом Python

Вы можете использовать это, сначала создав дерево k-d для ваших данных xyz:

import numpy as np  #for later code
from scipy.spatial import cKDTree

kdtree = cKDTree(xyzData)

Затем вы должны запросить дерево k-d с точкой point, чтобы вычислить расстояние между point и его ближайшим соседом. Результатом этого запроса является расстояние NN_dist между point и его ближайшим соседом и индекс NN_idx этого соседа. Чтобы вычислить это для всех ваших точек, нам нужен цикл for, но, учитывая алгоритм дерева k-d, это намного быстрее, чем вычисление методом грубой силы:

NN_dists = np.zeros(numPoints)  #pre-allocate an array to store distances
for i in range(numPoints):
    point = xyzData[i]

    NN_dist, NN_idx = kdtree.query(point,k=[1])

    #Note: 'k' specifies the kth neighbor distance to compute, 
    #so set k=2 if you end up finding the point as its own "neighbor":
    if NN_dist == 0:
        NN_dist, NN_idx = targetTree.query(curCoord,k=[2])

    NN_dists[i] = NN_dist

(подробнее см. запрос дерева k-d ).

Затем, чтобы найти расстояния, которые ниже некоторого порога, вы можете использовать встроенную утилиту массивов NumPy при использовании операторов сравнения (например, <):

distanceThres = 10
goodIdx = NN_dists < distanceThres
goodPoints = xyzData[goodIdx]

Это даст вам индексы goodIdx и точки goodPoints, которые находятся в пределах указанного вами порогового значения расстояния distanceThres (хотя вам, возможно, придется изменить этот код в зависимости от формы / формата ваших координатных данных xyz).

---

* Светлый фон на деревьях kd (приукрашивание мелких деталей - подробности см. В ссылках): метод дерева kd разбивает набор данных таким образом, чтобы избежать вычисления расстояния между каждой отдельной точкой (т. Е. Методом грубой силы) ). Это делается путем разделения набора данных на двоичные разделы пространства для создания дерева k-d. Эти разделы таковы, что вычисление расстояния (например, поиск ближайшего соседа) может игнорировать точки данных, которые находятся в удаленных разделах. Кроме того, это же дерево k-d повторно используется для каждой точки.

В Интернете много ресурсов о k-d деревьях в целом. Я нашел эти ссылки наиболее полезными, когда узнал об этом алгоритме: Стэнфордские k-d деревья или Принстонские k-d деревья .

Дайте мне знать, если у вас есть вопросы - у меня была именно эта проблема во время астрономического проекта, так что я могу помочь больше.

Answer 2

У меня нет прямого опыта с этим, но scipy.spatial.distance.pdist может быть тем, что вы ищете.

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