返回给定numpy数组的searchlight向量

Question

考虑一个维度为(30x40x50)的3Dnumpy数组D。对于每个体素D[x,y,z]，我希望存储一个包含某个半径内的相邻体素的向量(包括D[x,y,z]本身)。

(作为示例，这里是半径为2的球体的图片：https://puu.sh/wwIYW/e3bd63ceae.png)

有没有一种简单而快速的方法来编码呢？

我已经为它写了一个函数，但它非常慢，而且空闲最终会崩溃，因为我存储向量的数据结构变得太大了。

当前代码：

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

这里有一个方法可以做到这一点。为了提高效率，因此需要使用ndarray:这只考虑了完整的体素。边缘必须“手动”管理。

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)

所有代价高昂的操作都是矢量化的。在外部循环中，为了清晰起见，最好使用理解列表。

现在可以对体素执行矢量化操作。例如，最亮的体素：

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

当然，所有这些都是在内存中广泛存在的。您必须为大数据或体素划分空间。

Answer 2

由于您说数据结构太大，所以您可能必须为给定的体素动态计算向量。不过，您可以非常快速地完成此操作：

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]]

要在给定的数组上使用对象，您可以简单地执行以下操作：

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

这应该会给你相同的向量，你可以从

M_out[i,j,k]

在您的示例代码中，无需将所有结果一次存储在内存中。

这也可能进一步优化，特别是在坐标过滤方面，但这可能不是必需的。希望这能有所帮助。