Suchscheinwerfervektoren für ein gegebenes numpy Array zurückgeben

Question

Betrachten Sie ein 3D-Numy-Array D der Dimension, sagen wir, (30 x 40 x 50). Für jedes Voxel D[x,y,z] möchte ich einen Vektor speichern, der benachbarte Voxel innerhalb eines bestimmten Radius enthält (einschließlich der D[x,y,z] selbst).

(hier als Beispiel ist ein Bild einer solchen Kugel mit dem Radius 2: https://puu.sh/wwIYW/e3bd63ceae.png)

Gibt es eine einfache und schnelle Möglichkeit, dies zu codieren?

Ich habe eine Funktion dafür geschrieben, aber es ist schmerzhaft langsam und IDLE schließlich stürzt ab, weil die Datenstruktur, in der ich die Vektoren speichern soll, zu groß wird.

Aktuelle Code:

def searchlight(M_in): 
    radius = 4 
    [m,n,k] = M_in.shape 
    M_out = np.zeros([m,n,k],dtype=object) 
    count = 0 
    for i in range(m): 
     for j in range(n): 
      for z in range(k): 
       i_interval = list(range((i-4),(i+5))) 
       j_interval = list(range((j-4),(j+5))) 
       z_interval = list(range((z-4),(z+5)))     
       coordinates = list(itertools.product(i_interval,j_interval,z_interval)) 
       coordinates = [pair for pair in coordinates if ((abs(pair[0]-i)+abs(pair[1]-j)+abs(pair[2]-z))<=radius)] 
       coordinates = [pair for pair in coordinates if ((pair[0]>=0) and (pair[1]>=0) and pair[2]>=0) and (pair[0]<m) and (pair[1]<n) and (pair[2]<k)] 
       out = [] 
       for pair in coordinates: 
        out.append(M_in[pair[0],pair[1],pair[2]]) 
       M_out[i,j,z] = out 
       count = count +1 
    return M_out 

Answer 1

Hier eine Möglichkeit, das zu tun. Aus Gründen der Effizienz müssen Sie daher ndarrays verwenden: Hierbei werden nur vollständige Voxel berücksichtigt. Kanten müssen "von Hand" verwaltet werden.

from pylab import * 
a=rand(100,100,100) # the data 
r=4 
ra=range(-r,r+1) 

sphere=array([[x,y,z] for x in ra for y in ra for z in ra if np.abs((x,y,z)).sum()<=r]) 
# the unit "sphere" 

indcenters=array(meshgrid(*(range(r,n-r) for n in a.shape),indexing='ij')) 
# indexes of the centers of the voxels. edges are cut. 

all_inds=(indcenters[newaxis].T+sphere.T).T 
#all the indexes. 

voxels=np.stack([a[tuple(inds)] for inds in all_inds],-1) 
# the voxels. 
#voxels.shape is (92, 92, 92, 129) 

Alle kostspieligen Operationen sind vektorisiert. Verständnislisten sind für die Klarheit in der externen Schleife bevorzugt.

Sie können nun vektorisierte Operationen an Voxeln ausführen. Zum Beispiel das hellste Voxel:

light=voxels.sum(-1) 
print(np.unravel_index(light.argmax(),light.shape)) 
#(33,72,64) 

All dies ist natürlich umfangreich im Speicher. Sie müssen Ihren Platz für große Daten oder Voxel teilen.

Answer 2

Da Sie die Datenstruktur sagen zu groß ist, werden Sie wahrscheinlich den Vektor im Fluge für ein gegebenes Voxel zu berechnen haben. Sie können dies ziemlich schnell aber:

class SearchLight(object): 
    def __init__(self, M_in, radius): 
     self.M_in = M_in 
     m, n, k = self.M_in.shape 

     # compute the sphere coordinates centered at (0,0,0) 
     # just like in your sample code 
     i_interval = list(range(-radius,radius+1)) 
     j_interval = list(range(-radius,radius+1)) 
     z_interval = list(range(-radius,radius+1)) 
     coordinates = list(itertools.product(i_interval,j_interval,z_interval)) 
     coordinates = [pair for pair in coordinates if ((abs(pair[0])+abs(pair[1])+abs(pair[2]))<=radius)] 

     # store those indices as a template 
     self.sphere_indices = np.array(coordinates) 

    def get_vector(self, i, j, k): 

     # offset sphere coordinates by the requested centre. 
     coordinates = self.sphere_indices + [i,j,k] 
     # filter out of bounds coordinates 
     coordinates = coordinates[(coordinates >= 0).all(1)] 
     coordinates = coordinates[(coordinates < self.M_in.shape).all(1)] 

     # use those coordinates to index the initial array. 
     return self.M_in[coordinates[:,0], coordinates[:,1], coordinates[:,2]] 

Um das Objekt auf einem bestimmten Array zu verwenden, können Sie einfach tun:

sl = SearchLight(M_in, 4) 
# get vector of values for voxel i,j,k 
vector = sl.get_vector(i,j,k) 

Dies sollten Sie den gleichen Vektor geben Sie

M_out[i,j,k] 

bekommen würde

in Ihrem Beispielcode, ohne alle Ergebnisse auf einmal im Speicher zu speichern.

Dies kann wahrscheinlich auch weiter optimiert werden, insbesondere in Bezug auf die Koordinatenfilterung, aber es ist möglicherweise nicht notwendig. Ich hoffe, das hilft.