Способ вычисления точек для поверхности сферы быстрее, чем вложенная петля?

Question

Привет, ребята, на данный момент у меня есть массив numpy, представляющий Voxelgrid.Моя функция получает координаты 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 Methods использует маску, но это тоже медленно .. Есть ли лучший способ?и да, мой полярный угол должен идти только от 2 / 3pi - пи ..

Answer 1

В соответствии с python - Как мне отойти от школы мысли «для петли»?- Software Engineering Stack Exchange :

Обычно обработка данных в Python лучше всего выражается в виде итераторов ... Но NumPy выворачивает все наизнанку: лучший подход - это выразить алгоритмкак последовательность операций с целым массивом, чтобы минимизировать количество времени, затрачиваемое на медленный интерпретатор Python, и максимально увеличить количество времени, затрачиваемое на быстрые скомпилированные процедуры NumPy.

• Работа изнутри:начните с самого внутреннего цикла и посмотрите, можно ли векторизоватьзатем, когда вы это сделаете, выйдите на один уровень и продолжайте.
• (сохраните исходную версию функции (которая, как вы уверены, верна), чтобы вы могли проверить ее на соответствие улучшенным версиям как на корректностьи скорость.)

---

Итак, давайте посмотрим:

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 )

Каждое значение необходимо использовать повторно для каждого значения phi, поэтому numpy.repeat по длине PHI.

---

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

=>

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

необходимо повторять для каждого значения eta, поэтому numpy.tile по длине 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()