Wie Geschwindigkeit der Ablagerung Werte auf einem großen 3D-Gitter zu erhöhen

Question

Ich habe eine große Reihe von Datenpunkten in 3 Spalten Vektoren. Es gibt 10 Millionen Punkte mit x, y, z Koordinaten.

Ich voxelisiere diese Punkte (Zuordnung zu einem diskreten Gitter basierend auf der Belegung). Es gibt zwei Möglichkeiten, die Voxelisierung zu erreichen. Der erste Weg ist eine einfache Binning-Prozedur. Wenn der Punkt in eine bestimmte Bin fällt, steigt die Intensität der Bin um 1. Die andere Möglichkeit besteht darin, einen Punkt mehreren Bins zuzuordnen und die Intensität basierend auf der Entfernung von den Bin-Zentren zu erhöhen. Ich möchte die zweite Methode der Voxelisierung durchführen.

Ein einfaches 2D-Beispiel hierfür ist: Sagen Sie bitte Punkt x haben, y = 1.7,2.2 und haben einen gleichmäßig beabstandeten Raster mit dem Abstand zwischen den Knoten in 0,5 x und y.

Verfahren unter Verwendung von 1: Der Punkt würde x binned bekommen, y = 1.5,2 mit Intensität = 1

Verfahren unter Verwendung von 2: Der Punkt würde erhalten verteilt (x, y), (x- .5, y), (x + .5, y), (x, y - .5), (x, y + .5) Mit Intensitäten = (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] 

Dann danach plane ich durch diese riesige Array zurückzuzuverfolgen und all diese Intensitäten bis zu den gleichen Punkten in dem meine letzte Karte zu bekommen, aber das spielt keine Rolle, in dem mo ment.

Dieser Code funktioniert voxelize die 3D-Punkte, aber es ist super langsam. Gibt es eine Möglichkeit, dies weniger naiv und schneller zu machen? Ich dachte daran, vorher ein Gitter mit Koordinaten und Intensität von 0 an jedem Punkt zu erstellen und in einer Schleife die Intensität mit + = oder etwas ähnlichem zu aktualisieren.

Answer 1

Es ist möglich, for-Schleife zu entfernen und numpy Operationen zu verwenden, um es zu kümmern. Gleicher Code wie bei Ihnen, ohne Schleife und Indizierung ist ~ 60x schneller:

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 

Ich nehme an, die einige weitere Optimierung ist möglich, mit mehr Matrix orientierte Berechnung, aber ich habe es nicht fassen genau das, was Code zu tun hat: -/Wahrscheinlich werden die Grid-Werte [:, 1-7, 0-2] vor der Berechnung der Entfernung und dann, wenn nur die Grid-Werte für die Entfernungsberechnung verwendet werden, durch das Weglassen einiger nicht benötigter Zuweisungen etwas Zeit eingespart.