По заданным точкам сетки XYZ рассчитать объем сферы внутри

Question

У меня есть большая трехмерная сетка (~ 800 000 точек) равномерно расположенных xyz декартовых точек, и я хочу найти объем сферы внутри, основываясь на количестве точек, занимающих сферу. В настоящее время я использую scipy cKDTree и определяю все точки в пределах определенного радиуса начала координат (0,0,0) с помощью query_ball_point для оценки объема, но этот объем (grid_vol, ниже) часто сильно отличается (ошибка 50% или больше ) чем истинный объем (сфера_вол).

import numpy as np
from scipy import spatial
import math
#constants
xl = -3.15
xr =  1.75
yl = -2.0
yr =  2.0
zl = -1.15
zr = 3.9
spacing = 0.05
R = 3.5
cube = spacing ** 3  
#Create grid
x=np.arange(xl,xr,spacing)
y=np.arange(yl,yr,spacing)
z=np.arange(zl,zr,spacing)
x2,y2,z2=np.meshgrid(x,y,z,indexing='ij')
all_grid=np.array([x2.flatten(),y2.flatten(),z2.flatten()]).T
cube = spacing ** 3                  
point_tree = spatial.cKDTree(all_grid)                  # allgrid = evenly spaced rectangular grid
n_voxel = len(point_tree.query_ball_point((0,0,0), R))   # number of points in grid occupying sphere of radius R
grid_vol = n_voxel * cube                               # volume based on grid points
sphere_vol = 4 / 3 * math.pi * R ** 3                # vol of sphere w/ radius R to compare

в этом случае:

grid_vol = 78.1275
sphere_vol = 179.5944

Мне было интересно, существует ли известный модуль, который подойдет для этого применения измерения сферы из точек сетки

Answer 1

Эй, Гилиан, я попробовал твой код, и, насколько я могу судить, он работает просто отлично (некоторая ошибка из-за аппроксимации, но не большая

import numpy as np
from scipy import spatial
#define constants
l=-1
r=1
spacing=0.05
R=0.5
cube = spacing ** 3  
#Create grid
x=np.arange(l,r,spacing)
y=np.arange(l,r,spacing)
z=np.arange(l,r,spacing)
x2,y2,z2=np.meshgrid(x,y,z,indexing='ij')
all_grid=np.array([x2.flatten(),y2.flatten(),z2.flatten()]).T
# your code            
point_tree = spatial.cKDTree(all_grid)
n_voxel = len(point_tree.query_ball_point((0,0,0), R))   # number of points in grid occupying sphere of radius R
grid_vol = n_voxel * cube                               # volume based on grid points
sphere_volume = 4 / 3 * np.pi * R ** 3
print(grid_vol) # 0.5185000000000001
print(sphere_volume) # 0.5235987755982988