Как увеличить скорость внесения значений в большую 3D-сетку

У меня есть большой набор точек данных в 3 векторах-столбцах. Есть 10 миллионов точек с координатами x, y, z.

Я вокселизирую эти точки (назначая их дискретной сетке в зависимости от занятости). Существует два способа выполнить вокселизацию. Первый способ — это простая процедура бинирования, при которой, если точка попадает в определенный бин, интенсивность этого бина увеличивается на 1. Другой способ — назначить точку нескольким бинам и увеличить интенсивность в зависимости от расстояния от центров бинов. Я хочу выполнить второй метод вокселизации.

Простой двумерный пример этого: скажем, у вас есть точка x, y = 1,7, 2,2 и равномерная сетка с расстоянием 0,5 между узлами в x и y.

Используя метод 1: точка будет разделена на x, y = 1,5, 2 с интенсивностью = 1.

Используя метод 2: точка будет распределена на (x,y),(x-.5,y),(x+.5,y),(x,y-.5),(x,y+.5) с интенсивности=(distTOpoint1/sumDistances),(distTopoint2/sumDistances),...,(distTopoint5/sumDistances)

def floorPartial (value, resolution):
    return np.floor (value / resolution) * resolution 

def EucDistSq(x1,y1,z1,x2,y2,z2):
        return (x1-x2)**2+(y1-y2)**2+(z1-z2)**2

xCoord=100*np.random.random(10000000)
yCoordC=100*np.random.random(10000000)
zCoord=100*np.random.random(10000000)

Xspacing=.1
Yspacing=.1
zspacing=.1
Grid=np.empty([len(xCoord),8,4])

for i in range(len(xCoord)):

    Grid[i,0,:]=[xCoord[i],yCoordC[i],zCoord[i],0] #Save original Point

    #calculate voxel which it would go to if it was simple binning
    vX=floorPartial(xCoord[i],Xspacing) 
    vY=floorPartial(yCoordC[i],Yspacing)
    vZ=floorPartial(zCoord[i],Zspacing)


    d1=EucDistSq(xCoord[i],yCoordC[i],zCoord[i],vX,vY,vZ)
    d2=EucDistSq(xCoord[i],yCoordC[i],zCoord[i],vX+Xspacing,vY,vZ)
    d3=EucDistSq(xCoord[i],yCoordC[i],zCoord[i],vX-Xspacing,vY,vZ)
    d4=EucDistSq(xCoord[i],yCoordC[i],zCoord[i],vX,vY+Yspacing,vZ)
    d5=EucDistSq(xCoord[i],yCoordC[i],zCoord[i],vX,vY-Yspacing,vZ)
    d6=EucDistSq(xCoord[i],yCoordC[i],zCoord[i],vX,vY,vZ+Zspacing)
    d7=EucDistSq(xCoord[i],yCoordC[i],zCoord[i],vX,vY,vZ-Zspacing)

    dt=np.sum([d1,d2,d3,d4,d5,d6,d7])

    #VoxelX,VoxelY,VoxelZ,intensity
    Grid[i,1,:]=[vX,vY,vZ,d1/dt]
    Grid[i,2,:]=[vX+Xspacing,vY,vZ,d2/dt]
    Grid[i,3,:]=[vX-Xspacing,vY,vZ,d3/dt]
    Grid[i,4,:]=[vX,vY+Yspacing,vZ,d4/dt]
    Grid[i,5,:]=[vX,vY-Yspacing,vZ,d5/dt]
    Grid[i,6,:]=[vX,vY,vZ+Zspacing,d6/dt]
    Grid[i,7,:]=[vX,vY,vZ-Zspacing,d7/dt]

Затем, после этого, я планирую вернуться через этот гигантский массив и добавить все эти интенсивности для тех же точек, чтобы получить окончательную карту, но на данный момент это не имеет значения.

Этот код работает для вокселизации трехмерных точек, но он очень медленный. Есть ли способ сделать это менее наивно и быстрее? Я думал заранее создать сетку с координатами и интенсивностью 0 в каждой точке и в цикле обновить интенсивность с помощью += или что-то в этом роде.