Поиск индекса ближайшей точки в массивах numpy координат x и y
У меня есть два массива 2d numpy: x_array содержит позиционную информацию в направлении x, y_array содержит позиции в направлении y.
Затем у меня есть длинный список x, y точек.
Для каждой точки списка мне нужно найти индекс массива местоположения (указанного в массивах), который ближе всего к этой точке.
Я наивно создал некоторый код, который работает, основываясь на этом вопросе:
Найти ближайшее значение в массиве numpy
то есть.
import time
import numpy
def find_index_of_nearest_xy(y_array, x_array, y_point, x_point):
distance = (y_array-y_point)**2 + (x_array-x_point)**2
idy,idx = numpy.where(distance==distance.min())
return idy[0],idx[0]
def do_all(y_array, x_array, points):
store = []
for i in xrange(points.shape[1]):
store.append(find_index_of_nearest_xy(y_array,x_array,points[0,i],points[1,i]))
return store
# Create some dummy data
y_array = numpy.random.random(10000).reshape(100,100)
x_array = numpy.random.random(10000).reshape(100,100)
points = numpy.random.random(10000).reshape(2,5000)
# Time how long it takes to run
start = time.time()
results = do_all(y_array, x_array, points)
end = time.time()
print 'Completed in: ',end-start
Я делаю это над большим набором данных и очень хотел бы ускорить его.
Кто-нибудь может это оптимизировать?
Спасибо.
ОБНОВЛЕНИЕ: РЕШЕНИЕ после предложений @silvado и @justin (ниже)
# Shoe-horn existing data for entry into KDTree routines
combined_x_y_arrays = numpy.dstack([y_array.ravel(),x_array.ravel()])[0]
points_list = list(points.transpose())
def do_kdtree(combined_x_y_arrays,points):
mytree = scipy.spatial.cKDTree(combined_x_y_arrays)
dist, indexes = mytree.query(points)
return indexes
start = time.time()
results2 = do_kdtree(combined_x_y_arrays,points_list)
end = time.time()
print 'Completed in: ',end-start
Этот код ускорил мой код (поиск 5000 точек в матрицах 100x100) в 100 раз. Интересно, что использование scipy.spatial.KDTree(вместо scipy.spatial.cKDTree) дало сопоставимые тайминги моему наивному решению, поэтому определенно стоит использовать версию cKDTree...
Ответы
Ответ 1
scipy.spatial
также имеет реализацию дерева k-d: scipy.spatial.KDTree
.
Обычно подход заключается в том, чтобы сначала использовать точечные данные для создания дерева k-d. Вычислительная сложность этого порядка N log N, где N - количество точек данных. Запросы диапазона и поиск ближайших соседей могут быть выполнены с использованием сложности log N. Это намного эффективнее, чем просто циклирование по всем точкам (сложность N).
Таким образом, если у вас есть повторяющиеся запросы диапазона или ближайшего соседа, рекомендуется использовать дерево k-d.
Ответ 2
Вот пример scipy.spatial.KDTree
In [1]: from scipy import spatial
In [2]: import numpy as np
In [3]: A = np.random.random((10,2))*100
In [4]: A
Out[4]:
array([[ 68.83402637, 38.07632221],
[ 76.84704074, 24.9395109 ],
[ 16.26715795, 98.52763827],
[ 70.99411985, 67.31740151],
[ 71.72452181, 24.13516764],
[ 17.22707611, 20.65425362],
[ 43.85122458, 21.50624882],
[ 76.71987125, 44.95031274],
[ 63.77341073, 78.87417774],
[ 8.45828909, 30.18426696]])
In [5]: pt = [6, 30] # <-- the point to find
In [6]: A[spatial.KDTree(A).query(pt)[1]] # <-- the nearest point
Out[6]: array([ 8.45828909, 30.18426696])
#how it works!
In [7]: distance,index = spatial.KDTree(A).query(pt)
In [8]: distance # <-- The distances to the nearest neighbors
Out[8]: 2.4651855048258393
In [9]: index # <-- The locations of the neighbors
Out[9]: 9
#then
In [10]: A[index]
Out[10]: array([ 8.45828909, 30.18426696])
Ответ 3
Если вы можете массировать данные в нужном формате, быстрый способ - использовать методы в scipy.spatial.distance
:
http://docs.scipy.org/doc/scipy/reference/spatial.distance.html
В частности, pdist
и cdist
обеспечивают быстрые способы вычисления попарных расстояний.