这是一个使用 Matt Timmermans 示例的 python 示例。如果您有多个球体，则使用 numba jits 函数和并行 prange 会非常快。为任何起始坐标添加了体素。网格只是在 3d 网格中球体素所在的位置放置 1。还添加了一个可选的填充步骤。注意：这仅适用于非负坐标，因此对于使用体素坐标的所有体素，球体需要位于正空间中。
#/usr/bin/env python
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d
from numba.experimental import jitclass
from numba import types, njit, prange
spec = [
('radius', types.float64),
('svoxel', types.int64[:]),
('grid', types.b1[:,:,:])
]
@jitclass(spec)
class RasterizeSphere(object):
def __init__(self, svoxel, radius, grid):
self.grid = grid
self.svoxel = svoxel
self.radius = radius
R2 = np.floor(self.radius**2)
zx = np.int64(np.floor(self.radius))
x = 0
while True:
while x**2 + zx**2 > R2 and zx >= x:
zx -= 1
if zx < x:
break
z = zx
y = 0
while True:
while x**2 + y**2 + z**2 > R2 and z >= x and z >= y:
z -= 1
if z < x or z < y:
break
self.fill_all(x, y, z)
###### Optional fills the inside of sphere as well. #######
for nz in range(z):
self.fill_all(x, y, nz)
y += 1
x += 1
def fill_signs(self, x, y, z):
self.grid[x + self.svoxel[0], y + self.svoxel[1], z + self.svoxel[2]] = True
while True:
z = -z
if z >= 0:
y = -y
if y >= 0:
x = -x
if x >= 0:
break
self.grid[x + self.svoxel[0], y + self.svoxel[1], z + self.svoxel[2]] = True
def fill_all(self, x, y, z):
self.fill_signs(x, y, z)
if z > y:
self.fill_signs(x, z, y)
if z > x and z > y:
self.fill_signs(z, y, x)
@njit(parallel=True, cache=True)
def parallel_spheres(grid):
# prange just to show the insane speedup for large number of spheres disable visualization below if using large amount of prange.
for i in prange(2):
radius = 4.5
svoxel = np.array([5, 5, 5])
max = np.int64(np.ceil(radius**2))
rs = RasterizeSphere(svoxel, radius, grid)
points = np.where(rs.grid)
return np.array([*points])
def main():
# Make max large enough to hold the spheres.
max = 100
points = parallel_spheres(np.zeros((max, max, max), dtype=bool))
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.scatter3D(*points)
plt.show()
if __name__ == '__main__':
main()