如何提高将值存放到大型三维网格上的速度

Question

我在3个列向量中有一大组数据点。有1000万个点，有x，y，z坐标。

我正在对这些点进行体素化（根据占用情况将它们分配给离散网格）。有两种方法可以完成体素化。第一种方式是简单的分箱程序，如果该点落入某个箱内，则箱的强度增加1.另一种方法是将一个点分配给多个箱，并根据与箱中心的距离增加强度。我希望完成第二种体素化方法。

一个简单的2d示例是：假设您有点x，y = 1.7,2.2并且在x和y中的节点之间具有均匀间隔的网格，距离为.5。

使用方法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]

然后在此之后我计划回溯这个巨大的阵列并将所有这些强度添加到相同的点以获得我的最终地图，但这在目前无关紧要。

此代码用于对3d点进行体素化，但它非常慢。有没有办法不那么天真和快速地做到这一点？我想在每个点预先创建一个网格，坐标和强度为0，在循环中用+ =或类似的东西更新强度。

Answer 1

可以删除for循环并使用numpy操作来处理它。与没有for循环和索引的代码相同的代码快〜60倍：

def ver_2(xCoord, yCoord, zCoord, xSpacing, ySpacing, zSpacing):
    Grid = numpy.empty([len(xCoord), 8, 4])
    Grid[:, 0, 0] = xCoord
    Grid[:, 0, 1] = yCoord
    Grid[:, 0, 2] = zCoord
    Grid[:, 0, 3] = 0
    #
    vX = floorPartial(xCoord, xSpacing)
    vY = floorPartial(yCoord, ySpacing)
    vZ = floorPartial(zCoord, zSpacing)
    #
    d1 = EucDistSq(xCoord, yCoord, zCoord, vX, vY, vZ)
    d2 = EucDistSq(xCoord, yCoord, zCoord, vX+xSpacing, vY, vZ)
    d3 = EucDistSq(xCoord, yCoord, zCoord, vX-xSpacing, vY, vZ)
    d4 = EucDistSq(xCoord, yCoord, zCoord, vX, vY+ySpacing, vZ)
    d5 = EucDistSq(xCoord, yCoord, zCoord, vX, vY-ySpacing, vZ)
    d6 = EucDistSq(xCoord, yCoord, zCoord, vX, vY, vZ+zSpacing)
    d7 = EucDistSq(xCoord, yCoord, zCoord, vX, vY, vZ-zSpacing)
    #
    dt = numpy.sum([d1, d2, d3, d4, d5, d6, d7], axis=0)
    # VoxelX,VoxelY,VoxelZ,intensity
    Grid[:, 1] = numpy.stack((vX, vY, vZ, d1/dt), axis=-1)
    Grid[:, 2] = numpy.stack((vX+xSpacing, vY, vZ, d2/dt), axis=-1)
    Grid[:, 3] = numpy.stack((vX-xSpacing, vY, vZ, d3/dt), axis=-1)
    Grid[:, 4] = numpy.stack((vX, vY+ySpacing, vZ, d4/dt), axis=-1)
    Grid[:, 5] = numpy.stack((vX, vY-ySpacing, vZ, d5/dt), axis=-1)
    Grid[:, 6] = numpy.stack((vX, vY, vZ+zSpacing, d6/dt), axis=-1)
    Grid[:, 7] = numpy.stack((vX, vY, vZ-zSpacing, d7/dt), axis=-1)
    return Grid

我想更多的优化是可能的，更多的矩阵导向计算，但我没有准确掌握代码必须做什么： - /可能在计算距离之前设置网格的[：，1-7,0-2]值，并且而仅使用网格值进行距离计算可以减少一些时间，因为省略了一些不需要的分配。