一种比嵌套循环更快地计算球体表面点的方法？

Question

嘿，伙计们，现在我有一个表示体素栅格的numpy数组。我的函数得到坐标x，y，z，半径和一个值。我想把这些值加到坐标上，坐标是给定半径的球面的一部分。我尝试了一些方法，但都很慢：

def spheric Surface (x, y, z, r, value):
    phi = 0
    while phi <= (2*math.pi):
        eta = math.pi * 2 / 3
        while eta <= math.pi:
            xx = x + r * math.sin(eta) * math.cos(phi)
            yy = y + r * math.sin(eta) * math.sin(phi)
            zz = z + r * math.cos(eta)
            xx = int(xx*resoultion+0.5)
            yy = int(yy*resolution+0.5)
            zz = int(zz*resolution+0.5)
            voxelGrid[xx][yy][zz] += value

            eta += 1/10 * math.pi
        phi += 1/10 * math.pi

第一种方法使用球坐标，半径越大，eta+=必须越小。。这个方法很慢。。你知道吗

def sphericSurface(x, y, z, r, value):
tol = 0.6

grenz = math.pi * 2 / 3
mask = (np.logical_and(np.logical_and((sx[:, None, None] - x) ** 2 + (sy[None, :, None] - y) ** 2 + (sz[None, None, :] - z) ** 2 <= (r + tol)**2,
                                      (sx[:, None, None] - x) ** 2 + (sy[None, :, None] - y) ** 2 + (sz[None, None, :] - z) ** 2 >= (r - tol)**2),
                       (sz[None, None, :] - z) <= (r*math.cos(grenz))))
x, y, z = np.where(mask==True)
z *= 2
voxelGrid[x,y,z] += value

Secon方法使用掩码，但这也很慢。。有没有更好的办法？是的，我的极角应该只有2/3pi。。你知道吗

Answer 1

根据python - How do I move away from the “for-loop” school of thought? - Software Engineering Stack Exchange：

Normally, data processing in Python is best expressed in terms of iterators... But NumPy turns all that inside out: the best approach is to express the algorithm as a sequence of whole-array operations, to minimize the amount of time spent in the slow Python interpreter and maximize the amount of time spent in fast compiled NumPy routines.

• Work from the inside out: that is, start with the innermost loop and see if can be vectorized; then when you've done that, move out one level and continue.
• (Keep the original version of the function (which you are confident is correct) so that you can test it against your improved versions both for correctness and speed.)

---

那么，让我们看看：

xx = x + r * math.sin(eta) * math.cos(phi)
yy = y + r * math.sin(eta) * math.sin(phi)
zz = z + r * math.cos(eta)
xx = int(xx*resoultion+0.5)
yy = int(yy*resolution+0.5)
zz = int(zz*resolution+0.5)
voxelGrid[xx][yy][zz] += value

=>；（大写字母表示向量）

voxelGrid = ceiling (
            [ x + r * sin (ETA) * cos (PHI) ,
              y + r * sin (ETA) * sin (PHI) ,
              z + r * cos (ETA) ] * resolution )

---

eta = math.pi * 2 / 3                
while eta <= math.pi:
    <...>
    eta += 1/10 * math.pi

=>

ETA = range ( pi*2/3, pi, pi*1/10 )

对于phiso ^{}的每个值，每个值都需要重复使用PHI。你知道吗

---

phi = 0
while phi <= (2*math.pi):    
    <...>
    phi += 1/10 * math.pi

=>

PHI = range ( 0.0, 2*pi, pi*1/10 )

需要为etaso ^{}的每个值重复ETA。你知道吗

---

最终导致：

# numpy.linspace seems better for your task than arange
PHI = np.arange(0.0, 2*np.pi, np.pi*1/10)
ETA = np.arange(np.pi*2/3, np.pi, np.pi*1/10)

ETA, PHI = np.repeat(ETA, PHI.shape[0]), np.tile(PHI, ETA.shape[0])

XX = x + r * np.sin(ETA) * np.cos(PHI)
YY = y + r * np.sin(ETA) * np.sin(PHI)
ZZ = z + r * np.cos(ETA)

voxelGrid = np.vstack((XX,YY,ZZ))
voxelGrid = np.ceil(resolution * voxelGrid)

# plot it if you want
#import matplotlib.pyplot, mpl_toolkits.mplot3d
#fig=matplotlib.pyplot.figure()
#ax=fig.add_subplot(111,projection='3d')
#ax.scatter(*voxelGrid)
#fig.show()