Suppose a 2d array is given as:
arr = array([[1, 1, 1],
[4, 5, 8],
[2, 6, 9]])
if point=array([1,1])
is given then I want to calculate the euclidean distance from all indices of arr
to point (1,1). The result should be
array([[1.41 , 1. , 1.41],
[1. , 0. , 1. ],
[1.41 , 1. , 1.41]])
For loop is too slow to do these computations. Is there any faster method to achieve this using numpy or scipy?
Thanks!!!
Approach #1
You can use scipy.ndimage.morphology.distance_transform_edt
-
def distmat(a, index):
mask = np.ones(a.shape, dtype=bool)
mask[index[0],index[1]] = False
return distance_transform_edt(mask)
Approach #2
Another with NumPy-native tools -
def distmat_v2(a, index):
i,j = np.indices(a.shape, sparse=True)
return np.sqrt((i-index[0])**2 + (j-index[1])**2)
Sample run -
In [60]: a
Out[60]:
array([[1, 1, 1],
[4, 5, 8],
[2, 6, 9]])
In [61]: distmat(a, index=[1,1])
Out[61]:
array([[1.41421356, 1. , 1.41421356],
[1. , 0. , 1. ],
[1.41421356, 1. , 1.41421356]])
In [62]: distmat_v2(a, index=[1,1])
Out[62]:
array([[1.41421356, 1. , 1.41421356],
[1. , 0. , 1. ],
[1.41421356, 1. , 1.41421356]])
Other proposed solution(s) :
# https://stackoverflow.com/a/61629292/3293881 @Ehsan
def norm_method(arr, point):
point = np.asarray(point)
return np.linalg.norm(np.indices(arr.shape, sparse=True)-point)
Using benchit
package (few benchmarking tools packaged together; disclaimer: I am its author) to benchmark proposed solutions.
In [66]: import benchit
In [76]: funcs = [distmat, distmat_v2, norm_method]
In [77]: inputs = {n:(np.random.rand(n,n),[1,1]) for n in [3,10,50,100,500,1000,2000,5000]}
In [83]: T = benchit.timings(funcs, inputs, multivar=True, input_name='Length')
In [84]: In [33]: T.plot(logx=True, colormap='Dark2', savepath='plot.png')
So, distmat_v2
seems to be doing really well, We can further improve on it, by leveraging numexpr
.
Extend to array of indices
We could extend the listed solutions to cover for the generic/bigger case of list/array of indices w.r.t. whom we need to get euclidean distances at rest of the positions, like so -
def distmat_indices(a, indices):
indices = np.atleast_2d(indices)
mask = np.ones(a.shape, dtype=bool)
mask[indices[:,0],indices[:,1]] = False
return distance_transform_edt(mask)
def distmat_indices_v2(a, indices):
indices = np.atleast_2d(indices)
i,j = np.indices(a.shape, sparse=True)
return np.sqrt(((i-indices[:,0])[...,None])**2 + (j-indices[:,1,None])**2).min(1)
Sample run -
In [143]: a = np.random.rand(4,5)
In [144]: distmat_indices(a, indices=[[2,2],[0,3]])
Out[144]:
array([[2.82842712, 2. , 1. , 0. , 1. ],
[2.23606798, 1.41421356, 1. , 1. , 1.41421356],
[2. , 1. , 0. , 1. , 2. ],
[2.23606798, 1.41421356, 1. , 1.41421356, 2.23606798]])
On top of @Divakar's good solutions, if you are looking for something abstract, you can use:
np.linalg.norm(np.indices(arr.shape, sparse=True)-point)
Note that it works with numpy 1.17+ (argument sparse
is added on the versions 1.17+ of numpy). Upgrade your numpy and enjoy.
In case you have older than 1.17 version of numpy , you can add dimensions to your point
by using this:
np.linalg.norm(np.indices(arr.shape)-point[:,None,None], axis=0)
output for point=np.array([1,1])
and given array in question:
[[1.41421356 1. 1.41421356]
[1. 0. 1. ]
[1.41421356 1. 1.41421356]]
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With