如何优化基于球坐标绘制体素化球体的For循环代码？

randPoints是包含10000组(x,y,z)坐标的列表，xx、yy、zz是meshgrid的输出，形状为(201,201,201)。当前通过以下循环向201×201×201的体素立方体n中填充10000个体素化球体：

for j in range(len(randPoints)):
    #prepare coordinate x,y,z 
    xc_small = randPoints[j][0]
    yc_small = randPoints[j][1]
    zc_small = randPoints[j][2]

    #just unit conversion
    xcm = xc_small*1e-6
    ycm = yc_small*1e-6
    zcm = zc_small*1e-6

    #calculate f to make sphere of radius 0.5e-6
    f =  ((xx-xcm)**2+(yy-ycm)**2+(zz-zcm)**2 <= (0.5e-6)**2)
    
    #add 1 sphere into cube n
    n = n+f

该循环执行10000次，耗时约10000秒，亟需加速计算。n最终是填充了10000个体素化球体的201×201×201体素立方体。

变量初始化详情：

x = np.arange(-50, 50, 0.5)
y = np.arange(-50, 50, 0.5)
z = np.arange(0, 2*50, 0.5)
xx, yy, zz   = np.meshgrid(x, y, z)
n = np.ones((int(201),int(201),int(201)))

我认为f的计算过程是性能瓶颈，尝试过两种优化方式但都未能解决问题：

• 向量化一次性计算：内存占用过高导致会话崩溃，代码如下：

#prepare all (x,y,z) into a vector for one-time computation
xc_small = np.array(randPoints_small)[:, 0]
yc_small = np.array(randPoints_small)[:, 1]
zc_small = np.array(randPoints_small)[:, 2]

#just unit conversion
xcm = xc_small*1e-6
ycm = yc_small*1e-6
zcm = zc_small*1e-6

#one-time computaion of f
f = ((xx[:, :, :, np.newaxis]-xcm)**2 + (yy[:, :, :, np.newaxis]-ycm)**2 + (zz[:, :, :, np.newaxis]-zcm)**2 <= (0.5e-6)**2)
n = np.sum(f, axis=-1)

• 循环并行化：计算速度并未得到提升。

内容的提问来源于stack exchange，提问作者HammJ