时间戳

首页 > 解决方案 > 需要使用 Numpy 在较大的 3D 阵列中评估较小的 3D 阵列

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

问题描述

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

似乎 Numpy 中的许多拆分功能都不能很好地工作，除非较小的立方体可以均匀地分成较大的立方体。试图通过思考过程，我制作了一个 48x48x48 的立方体，它的顺序是从 1 到 110,592。然后我正在考虑使用以下代码将其重塑为 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]],

等等等等

标签： pythonarraysnumpy

解决方案

解决方案

在线有此解决方案的实时版本，您可以自己尝试

您的原始问题有一个简单的（某种）解决方案，即在 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 的立方体，但你想要的是一个由更小的立方体组成的立方体。一个可以通过改变其步幅来转换为另一个。对于 48x48x48 的 dtype 数组int64，步幅最初将设置为(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（意味着它们共享数据，因此不需要复制），形状为(16,16,16,3,3,3). 可以通过将索引传递给来访问第ijk3x3 立方体：arrarr2

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 数组的具体问题。Retriding 也可用于查找所有重叠的 3x3x3 立方体，如cubecube上面的函数中所示。

推荐阅读