Python 将 N 维体积划分为统一的子体积

Question

我有一个对象占据了一个底层的 N 维正方形网格(由一个 numpy 数组表示)，因此只有 25% 的网格点被占用。每个 1x1x1x... N 立方体(即超立方体)在此网格中包含相同数量的此对象(仅位于此超立方体的某些顶点)。我有一个所有占用的网格点的坐标数组。任务是循环遍历此数组并提取每个 1x1x1... 超立方体的占用坐标，并将它们存储在新数组中以供进一步处理。

这种情况最好用例子来解释。考虑已选择基础网格的 3-D 情况，以便 1<=i,j,k<=4给出 2d numpy 数组: A = [ [1 1 1] [1 2 1] [1 3 1] [1 4 1] [2 1 1] [2 2 1] [2 3 1] [2 4 1] [3 1 1] [3 2 1] [3 3 1] [3 4 1] [4 1 1] [4 2 1] [4 3 1] [4 4 1] [1 1 2] [1 2 2] [1 3 2] [1 4 2] [2 1 2] [2 2 2] [2 3 2] [2 4 2] [3 1 2] [3 2 2] [3 3 2] [3 4 2] [4 1 2] [4 2 2] [4 3 2] [4 4 2] [1 1 3] [1 2 3] [1 3 3] [1 4 3] [2 1 3] [2 2 3] [2 3 3] [2 4 3] [3 1 3] [3 2 3] [3 3 3] [3 4 3] [4 1 3] [4 2 3] [4 3 3] [4 4 3] [1 1 4] [1 2 4] [1 3 4] [1 4 4] [2 1 4] [2 2 4] [2 3 4] [2 4 4] [3 1 4] [3 2 4] [3 3 4] [3 4 4] [4 1 4] [4 2 4] [4 3 4] [4 4 4] ]

在这种情况下我需要处理的 2d numpy 数组的示例是: B = [[1,1,1],[1,2,1],[1,3,1],[1,4,1],[2,2,1],[2,3,1], [2,4,1],[3,2,1],[3,3,1],[3,4,1],[4,1,1],[4,3,1],[1 ,1,2],[1,4,2],[2,1,2],[2,2,2],[2,4,2],[3,1,2],[3,2 ,2],[3,4,2],[4,1,2],[4,2,2],[4,3,2],[4,4,2],[1,1,3 ],[1,4,3],[2,1,3],[2,2,3],[2,3,3],[2,4,3],[3,1,3], [3,2,3],[3,3,3],[4,1,3],[4,2,3],[4,3,3],[4,4,3],[1 ,2,4],[1,3,4],[1,4,4],[2,1,4],[2,2,4],[2,3,4],[3,1 ,4],[3,2,4],[3,3,4],[4,3,4],[4,4,4]]

B仅是占用的网格点的数组。骑车经过 B ，我想获取与底层网格的每个立方体相关的占用坐标。立方体边缘由 1<=i,j,k<=2 定义, ( 3<=i<=4 ) & ( 1<=j,k<=2 ), 等。 , 我想提取数组 [[1,1,1],[1,1,2],[1,2,1],[2,1,2],[2,2,1],[ 2,2,2]] 来自 1<=i,j,k<=2并将其存储在一个新数组中以供进一步处理。然后重复此操作，直到考虑到此示例中的所有 8 个立方体。请注意，网格始终选择为均匀的，以便始终可以进行这种类型的划分。虽然我可以选择均匀的网格，但我无法控制对象占用的位置。

Numpy 数组切片不起作用，因为一般情况下，B 中的网格点不连续。 .我尝试使用“for”循环的代码来为立方体的边缘施加范围，但它很快变得太复杂，而且看起来不是正确的方法。然后是“numpy.where”函数；但条件的复杂性使其相当棘手。 “numpy.extract”也出现了类似的挑战。

Answer 1

如果您将网格作为 ndim+1 维坐标数组，如下所示

a = np.stack(np.mgrid[1:5, 1:5], -1)

那么这只是一个明智的 reshape 和转置的问题:

import itertools as it

# dimensions for reshape:
# split all coordinate axes into groups of two, leave the
# "dimension subscript" axis alone     
chop = *it.chain.from_iterable((i//2, 2) for i in a.shape[:-1]),

# shuffle for transpose:
# regroup axes into coordinates then pairs then subscript
ax_shuf = *(1 ^ np.r_[1:4*a.ndim-3:2] % (2*a.ndim-1)), -1

# put it together and flatten, i.e. join coordinates and join pairs
a.reshape(*chop, -1).transpose(*ax_shuf).reshape(np.prod(chop[::2]), -1, a.shape[-1])

结果:

array([[[1, 1],
        [1, 2],
        [2, 1],
        [2, 2]],

       [[1, 3],
        [1, 4],
        [2, 3],
        [2, 4]],

       [[3, 1],
        [3, 2],
        [4, 1],
        [4, 2]],

       [[3, 3],
        [3, 4],
        [4, 3],
        [4, 4]]])

更新

如果 B 是完整网格的无序子集，您可以执行以下技巧:

1. 减去列最小值以确保奇偶校验从偶数开始
2. floor 除以 2，每个感兴趣的立方体的角将折叠成一个点
3. 现在在行上使用 np.unique 来获取组

注意:如果需要，您可以将下面代码中的 argsort 替换为 Most efficient way to sort an array into bins specified by an index array? 中的解决方案之一，从而加快速度。 .

B = np.array(B)
unq, idx, cts = np.unique((B-B.min(0))//2, axis=0, return_inverse=True, return_counts=True)
if (cts[0]==cts).all():
    result = B[idx.argsort()].reshape(len(unq), cts[0], -1)
else:
    result = np.split(B[idx.argsort()], cts[:-1].cumsum())

结果:

array([[[1, 1, 1],
        [1, 2, 1],
        [2, 2, 1],
        [2, 2, 2],
        [2, 1, 2],
        [1, 1, 2]],

       [[2, 2, 3],
        [2, 1, 3],
        [1, 1, 3],
        [2, 1, 4],
        [2, 2, 4],
        [1, 2, 4]],

       [[1, 4, 2],
        [2, 4, 2],
        [1, 3, 1],
        [1, 4, 1],
        [2, 3, 1],
        [2, 4, 1]],

       [[1, 4, 4],
        [1, 3, 4],
        [2, 4, 3],
        [2, 3, 3],
        [1, 4, 3],
        [2, 3, 4]],

       [[4, 1, 1],
        [4, 1, 2],
        [3, 2, 1],
        [3, 2, 2],
        [4, 2, 2],
        [3, 1, 2]],

       [[4, 2, 3],
        [3, 2, 4],
        [3, 1, 3],
        [3, 2, 3],
        [3, 1, 4],
        [4, 1, 3]],

       [[4, 4, 2],
        [4, 3, 2],
        [3, 4, 2],
        [4, 3, 1],
        [3, 4, 1],
        [3, 3, 1]],

       [[4, 3, 3],
        [3, 3, 3],
        [4, 3, 4],
        [3, 3, 4],
        [4, 4, 3],
        [4, 4, 4]]])