我有一个三维网格坐标
x = linspace(0, Lx, Nx)
y = linspace(0, Ly, Ny)
z = linspace(0, Lz, Nz)
和我需要个指数点(即X [I],Y [J],Z [K])内的位置(X0,Y0,Z0)的某些半径R。 N_i可能相当大。我可以做一个简单的循环来找到我所需要的。
and I need to index points (i.e. x[i],y[j],z[k]) within some radius R of a position (x0,y0,z0). N_i can be quite large. I can do a simple loop to find what I need
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的?
but this quite slow. Is there a more natural/ faster way to do this type of operation with numpy?
这个怎么样:
import scipy.spatial as sp
x = np.linspace(0, Lx, Nx)
y = np.linspace(0, Ly, Ny)
z = np.linspace(0, Lz, Nz)
#Manipulate x,y,z here to obtain the dimensions you are looking for
center=np.array([x0,y0,z0])
#First mask the obvious points- may actually slow down your calculation depending.
x=x[abs(x-x0)<cutoff]
y=y[abs(y-y0)<cutoff]
z=z[abs(z-z0)<cutoff]
#Generate grid of points
X,Y,Z=np.meshgrid(x,y,z)
data=np.vstack((X.ravel(),Y.ravel(),Z.ravel())).T
distance=sp.distance.cdist(data,center.reshape(1,-1)).ravel()
points_in_sphere=data[distance<cutoff]
相反,最后两行,你应该能够做到:
Instead of the last two lines you should be able to do:
tree=sp.cKDTree(data)
mask=tree.query_ball_point(center,cutoff)
points_in_sphere=data[mask]
如果你不想叫空间:
distance=np.power(np.sum(np.power(data-center,2),axis=1),.5)
points_in_sphere=data[distance<cutoff]
上一篇:是否的P / Invoke在64位Windows需要不同的签名比对32位?比对、不同、Invoke、Windows
下一篇:GLSurfaceView onDrawFrame清算行为行为、GLSurfaceView、onDrawFrame