How to pair (x,y) pairs using numpy

Question

I have 2, 2D NumPy arrays consisting of ~300,000 (x,y) pairs each. Given the i^th (x, y) pair in array A, I need to find the corresponding j^th (x,y) pair in array B such that ([x_i - x_j]² + [y_i - y_j]²)^½, the distance between the two (x,y) pairs, is minimized.

What I'm currently doing is argmin like so:

thickness = []
for i in range(len(A)):
    xi = A[i][0]
    yi = A[i][1]
    idx = (np.sqrt(np.power(B[:, 0] - xi, 2) + np.power(B[:, 1] - yi, 2))).argmin()
    thickness.append([(xi + B[idx][0]) / 2, (yi + B[idx][1]) / 2,
                      A[i][2] + B[idx][2]])

Is there a faster way to accomplish this?

Answer 1

We can use Cython-powered kd-tree for quick nearest-neighbor lookup to get nearest neighbour and hence achieve our desired output, like so -

from scipy.spatial import cKDTree

idx = cKDTree(B[:,:2]).query(A[:,:2], k=1)[1]
thickness = [(A[:,0] + B[idx,0]) / 2, (A[:,1] + B[idx,1]) / 2, A[:,2] + B[idx,2]]

Answer 2

you can use numpy broadcast to find the indexes in B with minimum distance like below. With broadcast distance is calculated between each row of A and B

A = np.random.rand(100,3)
B = np.random.rand(90,3)

diff = A[:, np.newaxis, :2] - B[:,:2]
distance = np.sum(diff**2, axis=2)
indx = np.argmin(distance, axis=1)
result = (A+B[indx])/2
result

Answer 3

The fastest solution will require you to sort array B by their distance to the origin. You will need to break up array B into 4 hastables for each quadrant. You will need to find the A[Xi-Nx] [Yi-Ny] the index which is closest in distance to to the origin to our point B(X,Y) for each quadrant.Most of the time the closest point in A to our point in B would reside in the same quadrant. For the case in which the closest point is a different quadrant you will need to iterate through point in order of distance to the current quadrant. if you only need to find the distance to one point you could sort array B based on that distance to Initial point in A from smallest to largest. This would mean that the minimun distance will always be at index 0. THis solution will run in N(LOG(N)) complexity which is basically the cost of sorting.

import numpy as np
def quadrant(point): 
    x= point[0]
    y = point [1]
    if (x > 0 and y > 0): 
        return 1

    elif (x < 0 and y > 0): 
        return 2

    elif (x < 0 and y < 0): 
        return 3

    elif (x > 0 and y < 0): 
       return 4



    else: 
       return 0


def calculateDistance(point1,point2):
    #calculates the distance sqr
    x1 = point1[0]  
    y1 = point1[1]
    x2 = point2[0]  
    y2 = point2[1]
    dist = ((x2 - x1)**2 + (y2 - y1)**2)  
    return dist  

A = np.random.rand(10,2)
B = np.random.rand(10,2)

#sort b By distance to the origin and put it in sub_Quadran array -- 1 array per quadran
origin = [0,0]     
q1= {}  
q2= {}  
q3= {}  
q4= {}    


def quadrant(point): 
    x= point[0]
    y = point [1]
    if (x > 0 and y > 0): 
       q1 [calculateDistance(origin,point)] = point


    elif (x < 0 and y > 0): 
        q2 [calculateDistance(origin,point)] = point

    elif (x < 0 and y < 0): 
      q3  [calculateDistance(origin,point)] = point

    elif (x > 0 and y < 0): 
       q4  [calculateDistance(origin,point)] = point



    else: 
       print ("Point is the origin")




def find_near(point):
#for now assume point is on first quadrant
  dist = calculateDistance(point,origin) # get  to origin distance sqr
  # =find the point with closses distance this can be more efficient 

  idx=  min(q1_keys, key=lambda x:abs(x-dist))
  return q1[idx]




for  point in B:
    quadrant(point)


q1_keys = list()
q2_keys = list()
q3_keys = list()
q4_keys = list()
for i in q1.keys():
    q1_keys.append(i)
#Sorting can be done in the step bellow but done here to speedpup implementation
q1_keys.sort()

for i in q2.keys():
    q2_keys.append(i)
#Sorting can be done in the step bellow but done here to speedpup implementation
q2_keys.sort()

for i in q3.keys():
    q3_keys.append(i)
#Sorting can be done in the step bellow but done here to speedpup implementation
q3_keys.sort()

for i in q4.keys():
    q4_keys.append(i)
#Sorting can be done in the step bellow but done here to speedpup implementation
q4_keys.sort()



b_1 = [0.01, 0.01]   

q1_keys.sort() #point sorted as closest to the origin
print (q1)   
print(q1_keys)   
print (find_near(b_1))

` Another solution is to break the array B into multiple sub arrays and find the distance of the input point in A in those B sub_arrays. See algorithm details in the link bellow: https://courses.cs.washington.edu/courses/cse421/11su/slides/05dc.pdf