Numpy 将较小的矩阵添加到较大的矩阵

Question

我有大的 3D 矩阵，指示代理在 3D 空间中的位置。如果上面没有代理，则矩阵的值为 0，如果上面有代理，则为 1。然后，我的问题是我希望代理“成长”，因为我希望它们由一个立方体（3x3x3）决定。如果已经有办法做到这一点，但是当代理靠近边界时我遇到了麻烦。例如，我有一个位置矩阵 100x100x100，如果我知道我的代理在位置 (x, y, z) 我会这样做：

positions_matrix = numpy.zeros((100, 100, 100))
positions_matrix[x - 1: x + 2, y - 1: y + 2, z - 1: z + 2] +=  numpy.ones((3, 3, 3))

当然，在我的真实代码中，我会循环更多位置，但基本上就是这样。这是可行的，但是当代理接近无法求和的边界时，问题就来了，因为切片得到的矩阵将小于那些矩阵。

知道如何解决它，或者 numpy 或任何其他包是否有这个实现？虽然我很确定我不是第一个面对这个问题的人，但我找不到它。

Answer 1

解决问题的更程序化的方法：import numpy as npm = np.zeros((100, 100, 100))slicing = tuple(    slice(max(0, x_i - 1), min(x_i + 2, d - 1))    for x_i, d in zip((x, y, z), m.shape))ones_shape = tuple(s.stop - s.start for s in slicing)m[slicing] += np.ones(ones_shape)但在其他方面与接受的答案相同。

Answer 2

您应该使用以下方法在下限和上限处进行切割：import numpy as npm = np.zeros((100, 100, 100))x_min, x_max = np.max([0, x-1]), np.min([x+2, m.shape[0]-1])y_min, y_max = np.max([0, y-1]), np.min([y+2, m.shape[1]-1])z_min, z_max = np.max([0, z-1]), np.min([z+2, m.shape[2]-1])m[x_min:x_max, y_min:y_max, z_min:z_max] += np.ones((x_max-x_min, y_max-y_min, z_max-z_min))

Answer 3

有一个使用np.put, 及其'clip'选项的解决方案。它只需要一点技巧，因为该函数需要扁平矩阵中的索引；幸运的是，该功能np.ravel_multi_index完成了这项工作：import itertoolsimport numpy as npx, y, z = 2, 0, 4positions_matrix = np.zeros((100,100,100))indices = np.array( list( itertools.product( (x-1, x, x+1), (y-1, y, y+1), (z-1, z, z+1)) ))flat_indices = np.ravel_multi_index(indices.T, positions_matrix.shape, mode='clip')positions_matrix.put(flat_indices, 1+positions_matrix.take(flat_indices))# positions_matrix[2,1,4] is now 1.0这个解决方案的好处是您可以使用其他模式，例如'wrap'（如果您的代理住在甜甜圈上;-) 或周期性空间）。我将解释它如何在较小的 2D 矩阵上工作：import itertoolsimport numpy as nppositions_matrix = np.zeros((8,8))ones = np.ones((3,3))x, y = 0, 4indices = np.array( list( itertools.product( (x-1, x, x+1), (y-1, y, y+1) )))# array([[-1,  3],#        [-1,  4],#        [-1,  5],#        [ 0,  3],#        [ 0,  4],#        [ 0,  5],#        [ 1,  3],#        [ 1,  4],#        [ 1,  5]])flat_indices = np.ravel_multi_index(indices.T, positions_matrix.shape, mode='clip')# array([ 3,  4,  5,  3,  4,  5, 11, 12, 13])positions_matrix.put(flat_indices, ones, mode='clip')# positions_matrix is now:# array([[0., 0., 0., 1., 1., 1., 0., 0.],#        [0., 0., 0., 1., 1., 1., 0., 0.],#        [0., 0., 0., 0., 0., 0., 0., 0.],#        [ ...顺便说一句，在这种情况下mode='clip'是多余的put。好吧，我只是作弊put做作业。+=1需要take和put：_positions_matrix.put(flat_indices, ones.flat + positions_matrix.take(flat_indices))# notice that ones has to be flattened, or alternatively the result of take could be reshaped (3,3)# positions_matrix is now: # array([[0., 0., 0., 2., 2., 2., 0., 0.],#        [0., 0., 0., 2., 2., 2., 0., 0.],#        [0., 0., 0., 0., 0., 0., 0., 0.],#        [ ...与其他解决方案相比，此解决方案有一个重要区别：ones矩阵始终为 (3,3)，这可能是优势，也可能不是优势。诀窍就在这个 flat_indices 列表中，它有重复的条目（剪辑的结果）。因此，如果您在最大索引处添加非常数子矩阵，则可能需要采取一些预防措施：x, y = 1, 7values = 1 + np.arange(9)indices = np.array( list( itertools.product( (x-1, x, x+1), (y-1, y, y+1) )))flat_indices = np.ravel_multi_index(indices.T, positions_matrix.shape, mode='clip')positions_matrix.put(flat_indices, values, mode='clip')# positions_matrix is now:# array([[0., 0., 0., 2., 2., 2., 1., 3.],#        [0., 0., 0., 2., 2., 2., 4., 6.],#        [0., 0., 0., 0., 0., 0., 7., 9.],...您可能期望最后一列是 2 5 8。目前，您可以处理flat_indices，例如通过放置-1越界位置。但是，如果np.put接受非平面索引，或者如果有一个 clip ，这一切都会更容易mode='ignore'。