индексирование сферического подмножества данных трехмерной сетки в numpy
У меня есть 3d сетка с координатами
x = linspace(0, Lx, Nx)
y = linspace(0, Ly, Ny)
z = linspace(0, Lz, Nz)
и мне нужно индексировать точки (то есть x [i], y [j], z [k]) в пределах некоторого радиуса R позиции (x0, y0, z0). N_i может быть довольно большим. Я могу сделать простой цикл, чтобы найти то, что мне нужно
points=[]
i0,j0,k0 = floor( (x0,y0,z0)/grid_spacing )
Nr = (i0,j0,k0)/grid_spacing + 2
for i in range(i0-Nr, i0+Nr):
for j in range(j0-Nr, j0+Nr):
for k in range(k0-Nr, k0+Nr):
if norm(array([i,j,k])*grid_spacing - (x0,y0,k0)) < cutoff:
points.append((i,j,k))
но это довольно медленно. Есть ли более естественный / более быстрый способ сделать этот тип операции с NumPy?