大規模な3Dグリッドに値を入金する速度を上げる方法

Question

3つの列ベクトルに大きなデータセットがあります。 x、y、z座標には1000万点があります。

私はこれらの点をボクセル化（占有に基づいて離散グリッドに割り当てる）しています。ボクセル化を達成するには2つの方法があります。第1の方法は、ビンの強度が1だけ増加する特定のビン内にポイントがある場合の単純なビニング手順である。もう1つの方法は、複数のビンにポイントを割り当て、ビンセンターからの距離に基づいて強度を増加させることである。私はボクセル化の第2の方法を達成したい。

単純な2番目の例は次のとおりです。 ポイントx、y = 1.7,2.2があるとします。 xとyのノード間で距離が.5の等間隔のグリッドがあります。点が（X、（x、y）に分配さになるだろう ：方法2を使用し 点はxにビニングになるだろう、Yは強度で1.5,2を= = 1

：方法1を使用し

（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操作を使用して世話をすることは可能です。ループとインデックス作成のためのないあなたと同じコードさ〜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]の値を距離計算の前に設定し、グリッド値のみを距離計算に使用すると、不要な割り当てを省略するため時間が短縮される可能性があります。