np.where（）以消除坐标太近的数据

Question

我正在对星团进行孔径光度测定，并且为了更容易地检测背景信号，我只希望看距n个像素远的恒星（在我的情况下为n = 16）。我有2个数组xs和ys，所有星星坐标的x和y值：我应该使用np.where来找到所有恒星的索引，其中到所有其他恒星的距离> = n

到目前为止，我的方法一直是for-loop

import numpy as np

# Lists of coordinates w. values between 0 and 2000 for 5000 stars

xs = np.random.rand(5000)*2000
ys = np.random.rand(5000)*2000

# for-loop, wherein the np.where statement in  question is situated
n = 16
for i in range(len(xs)):
    index = np.where( np.sqrt( pow(xs[i] - xs,2) + pow(ys[i] - ys,2)) >= n)

由于星星密集地聚集在一起，所以我期望数据会大大减少，尽管即使当我尝试将n = 1000时，我仍然剩下大约4000个数据点

Answer 1

仅使用numpy（和部分答案here）

X = np.random.rand(5000,2) * 2000
XX = np.einsum('ij, ij ->i', X, X)
D_squared = (XX[:, None] + XX - 2 * X.dot(X.T)
out = np.where(D_squared.min(axis = 0) > n**2)

使用scipy.spatial.pdist

from scipy.spatial import pdist, squareform
D_squared = squareform(pdist(x, metric = 'sqeuclidean'))
out = np.where(D_squared.min(axis = 0) > n**2)

使用KDTree获得最大的速度：

from scipy.spatial import KDTree

X_tree = KDTree(X)
in_radius = np.array(list(X_tree.query_pairs(n))).flatten()
out = np.where(~np.in1d(np.arange(X.shape[0]), in_radius))

Answer 2

np.random.seed(seed=1)
xs = np.random.rand(5000,1)*2000
ys = np.random.rand(5000,1)*2000

n = 16

mask = (xs>=0)

for i in range(len(xs)):
    if mask[i]:
        index = np.where( np.sqrt( pow(xs[i] - x,2) + pow(ys[i] - y,2)) <= n)
        mask[index] = False
        mask[i] = True

x = xs[mask]
y = ys[mask]

print(len(x))

4220