CODEWITHSUNDEEP
pythonarraysnumpyscipydistance
fabian
fabian

Reputation: 1881

Find pair of index values that minimizes the Euclidian distance between two meshgrids and a column vector

I want to find the two find the pair of values and their index number in a meshgrid that a closets to another pair of values. Suppose I have two vectors a= np.array([0.01,0.5,0.9]) and b = np.array([0,3,6,10]) and two meshgrids X,Y = np.meshgrid(a,b). For illustration, they look as follows:

X= array([[ 0.1,  0.5,  0.9],
         [ 0.1,  0.5,  0.9],
         [ 0.1,  0.5,  0.9],
         [ 0.1,  0.5,  0.9]])

Y =array([[ 0,  0,  0],
         [ 3,  3,  3],
         [ 6,  6,  6],
         [10, 10, 10]])

Now, I have another array called c of dimension (2 x N). For illustration suppose c contains the following entries: 

c = array([[ 0.07268017,  0.08816632,  0.11084398,  0.13352165,  0.1490078 ],
           [ 0.00091219,  0.00091219,  0.00091219,  0.00091219,  0.00091219]])

Denote a column vector of c by x. For each vector x I want to find enter image description here

To complicate matters a bit, I am in fact not only looking for the index with the smallest distance (i,j) but also the second smallest distance (i',j').

All my approaches so far turned out to be extremely complicated and involved a lot of side routes. Does someone have an idea for how to tackle the problem efficiently?

Upvotes: 0

Views: 689

Answers (2)

gboffi
gboffi

Reputation: 25093

This is more a comment than an answer but i like to [... lots of stuff mercifully deleted, that you can still see using the revision history ...]

Complete Revision

As a followup of my own comment, please look at the following

Setup

In [25]: from numpy import *
In [26]: from scipy.spatial import KDTree
In [27]: X= array([[ 0.1,  0.5,  0.9],                                                
         [ 0.1,  0.5,  0.9],                                                  
         [ 0.1,  0.5,  0.9],
         [ 0.1,  0.5,  0.9]])
In [28]: Y =array([[ 0,  0,  0],                                                      
         [ 3,  3,  3],                                                        
         [ 6,  6,  6],
         [10, 10, 10]])
In [29]: c = array([[ 0.07268017,  0.08816632,  0.11084398,  0.13352165,  0.1490078 ],
           [ 0.00091219,  0.00091219,  0.00091219,  0.00091219,  0.00091219]])

Solution

Two lines of code, please notice that you have to pass the transpose of your c array.

In [30]: tree = KDTree(zip(X.ravel(), Y.ravel()))
In [31]: tree.query(c.T,k=2)                                           
Out[31]: 
(array([[ 0.02733505,  0.4273208 ],
        [ 0.01186879,  0.41183469],
        [ 0.01088228,  0.38915709],
        [ 0.03353406,  0.36647949],
        [ 0.04901629,  0.35099339]]), array([[0, 1],
        [0, 1],
        [0, 1],
        [0, 1],
        [0, 1]]))

Comment

To interpret the result, the excellent scipy docs inform you that tree.query() gives you back two arrays, containing respectively for each point in c

  1. a scalar or an array of length k>=2 giving you the distances from the point to the closest point on grid, the second closest, etc,
  2. a scalar or an array of length k>=2 giving you the indices pointing to the grid point(s) closest (next close etc).

To access the grid point, KDTree maintains a copy of the grid data, e.g

In [32]: tree.data[[0,1]]
Out[32]: 
array([[ 0.1,  0. ],
       [ 0.5,  0. ]])

where [0,1] is the first element of the second output array.

Should you need the indices of the closest(s) point in the mesh matrices, it simply a matter of using divmod.

Upvotes: 1

Frank M
Frank M

Reputation: 1580

If X, Y always come from meshgrid(), your minimization is separable in X and Y. Just find the closest elements of X to c[0,] and the closest elements of Y to c[1,] --- you don't need to calculate the 2-dimensional metric.

If either a or b have uniform steps, you can save yourself even more time if you scale the corresponding values of c onto the indexes. In your example, all(a == 0.1+0.4*arange(3)), so you can find the x values by inverting: x = (c[0,] - 0.1)/0.4. If you have an invertible (possibly non-linear) function that maps integers onto b, you can similarly find y values directly by applying the inverse function to c[1,].

Upvotes: 2

Related Questions