Znajdowanie indeksu najbliższego punktu w numpy tablicach współrzędnych x i y

Mam dwie tablice numpy 2d: x_array zawiera informacje o położeniu w kierunku x, y_array zawiera pozycje w kierunku y.

Mam wtedy długą listę punktów x, y.

Dla każdego punktu na liście muszę znaleźć indeks tablicy lokalizacji (określonej w tablicach), która jest najbliższa temu punktowi.

Naiwnie stworzyłem kod, który działa, w oparciu o następujące pytanie: Znajdź najbliższą wartość w tablicy numpy

to znaczy

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

Robię to na dużym zbiorze danych i naprawdę chciałbym to trochę przyspieszyć. Czy ktoś może to zoptymalizować?

Dzięki.

AKTUALIZACJA: ROZWIĄZANIE po sugestiach @silvado i @justin (poniżej)

# 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

Powyższy kod przyspieszył mój kod (wyszukiwanie 5000 punktów w macierzach 100x100) 100 razy. Co ciekawe, użycie scipy.spatial.KDTree (zamiast scipy.spatial.cKDTree ) dało porównywalne czasy do mojego naiwnego rozwiązania, więc zdecydowanie warto skorzystać z wersji cKDTree ...

scipy.spatial ma również implementację drzewa kd: scipy.spatial.KDTree .

Podejście polega na tym, że najpierw używa się danych punktów do zbudowania drzewa kd. Złożoność obliczeniowa tego jest rzędu N log N, gdzie N to liczba punktów danych. Zapytania o zasięg i wyszukiwanie najbliższych sąsiadów można następnie przeprowadzać ze złożonością dziennika N. Jest to o wiele bardziej wydajne niż zwykła jazda na rowerze przez wszystkie punkty (złożoność N).

Tak więc, jeśli masz powtarzające się zapytania dotyczące zasięgu lub najbliższego sąsiada, wysoce zalecane jest drzewo kd.

Oto scipy.spatial.KDTreeprzykład

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])

Jeśli możesz wmasować swoje dane w odpowiedni format, najszybszym sposobem jest użycie metod w scipy.spatial.distance:

http://docs.scipy.org/doc/scipy/reference/spatial.distance.html

W szczególności pdist i cdistzapewniają szybkie sposoby obliczania odległości parami.