python - 需要使用 Numpy 评估较大 3D 阵列中的较小 3D 阵列

Question

我必须取一个 50x50x50 的随机整数数组，并确定其中哪个连续的 3x3x3 立方体的总和最大。

Numpy 中的许多拆分功能似乎都无法正常工作，除非较小的立方体可以均匀地分成较大的立方体。尝试通过思考过程，我制作了一个 48x48x48 的立方体，它的顺序正好是从 1 到 110,592。然后我正在考虑使用以下代码将其 reshape 为 4D 数组并评估哪个数组的总和最大？当我输入此代码时，它会按照不理想的顺序拆分数组。我希望第一个数组是 3x3x3 立方体，该立方体位于 48x48x48 立方体的一角。我可以添加一种语法来实现这一点吗？

import numpy as np

arr1 = np.arange(0,110592)
arr2=np.reshape(arr1, (48,48,48))
arr3 = np.reshape(arr2, (4096, 3,3,3))
arr3

输出:

array([[[[     0,      1,      2],
         [     3,      4,      5],
         [     6,      7,      8]],

        [[     9,     10,     11],
         [    12,     13,     14],
         [    15,     16,     17]],

        [[    18,     19,     20],
         [    21,     22,     23],
         [    24,     25,     26]]],

期望的输出:

array([[[[     0,      1,      2],
         [    48,      49,      50],
         [     96,      97,      98]],

等等等等

Answer 1

解决方案

There's a live version of this solution online you can try for yourself

对于在 50x50x50 立方体中找到最大 3x3x3 子立方体的原始问题，有一个简单的(某种)解决方案，该解决方案基于更改输入数组的步幅。这个解决方案是完全矢量化的(意味着没有循环)，因此应该从 Numpy 中获得最佳性能:

import numpy as np

def cubecube(arr, cshape):
    strides = (*arr.strides, *arr.strides)
    shape = (*np.array(arr.shape) - cshape + 1, *cshape)
    return np.lib.stride_tricks.as_strided(arr, shape=shape, strides=strides)

def maxcube(arr, cshape):
    cc = cubecube(arr, cshape)
    ccsums = cc.sum(axis=tuple(range(-arr.ndim, 0)))
    ix = np.unravel_index(np.argmax(ccsums), ccsums.shape)[:arr.ndim]
    return ix, cc[ix]

maxcube 函数接受一个数组和子立方体的形状，并返回一个元组 (first-index-of-largest-cube, largest-cube) .下面是如何使用 maxcube 的示例:

shape = (50, 50, 50)
cshape = (3, 3, 3)

# set up a 50x50x50 array
arr = np.arange(np.prod(shape)).reshape(*shape)

# set one of the subcubes as the largest
arr[37, 26, 11] = 999999

ix, cube = maxcube(arr, cshape)
print('first index of largest cube: {}'.format(ix))
print('largest cube:\n{}'.format(cube))

哪些输出:

first index of largest cube: (37, 26, 11)
largest cube:
[[[999999  93812  93813]
  [ 93861  93862  93863]
  [ 93911  93912  93913]]

 [[ 96311  96312  96313]
  [ 96361  96362  96363]
  [ 96411  96412  96413]]

 [[ 98811  98812  98813]
  [ 98861  98862  98863]
  [ 98911  98912  98913]]]

深入解释

立方体的立方体

您拥有的是一个 48x48x48 的立方体，但您想要的是由更小的立方体组成的立方体。一个可以通过改变它的步幅转换成另一个。对于 dtype int64 的 48x48x48 数组，步幅最初将设置为 (48*48*8, 48*8, 8)。每个非重叠的 3x3x3 子立方体的第一个值可以用 (3*48*48*8, 3*48*8, 3*8) 步幅迭代。结合这些步幅来获得立方体立方体的步幅:

# Set up a 48x48x48 array, like in OP's example
arr = np.arange(48**3).reshape(48,48,48)

shape = (16,16,16,3,3,3)
strides = (3*48*48*8, 3*48*8, 3*8, 48*48*8, 48*8, 8)

# restride into a 16x16x16 array of 3x3x3 cubes
arr2 = np.lib.stride_tricks.as_strided(arr, shape=shape, strides=strides)

arr2 是 arr 的一个 View (意味着它们共享数据，所以不需要复制)，形状为 (16,16, 16,3,3,3)。可以通过将索引传递给 arr2 来访问 arr 中的第 ijk 个 3x3 立方体:

i,j,k = 0,0,0
print(arr2[i,j,k])

输出:

[[[   0    1    2]
  [  48   49   50]
  [  96   97   98]]

 [[2304 2305 2306]
  [2352 2353 2354]
  [2400 2401 2402]]

 [[4608 4609 4610]
  [4656 4657 4658]
  [4704 4705 4706]]]

您只需对内轴求和即可获得所有子立方体的总和:

sumOfSubcubes = arr2.sum(3,4,5)

这将产生一个 16x16x16 数组，其中每个值都是原始数组中一个非重叠的 3x3x3 子立方体的总和。这解决了 OP 询问的有关 48x48x48 数组的具体问题。 Restriding 也可用于查找所有重叠的 3x3x3 立方体，如上面的 cubecube 函数所示。