将 Python MeshGrid 拆分为单元格

Posted 2023-02-22

【中文标题】将 Python MeshGrid 拆分为单元格

【英文标题】：Splitting Python MeshGrid Into Cells

【发布时间】：2018-12-09 06:09:57

【问题描述】：

问题陈述

需要将 N 维 MeshGrid 拆分为“立方体”：

例如） 二维案例：

(-1,1) |(0,1) |(1,1)

(-1,0) |(0,0) |(1,0)

(-1,-1)|(0,-1)|(1,-1)

将有 4 个单元格，每个单元格有 2^D 个点：

我希望能够处理网格，将每个单元格的坐标点放入容器中进行进一步处理。

Cells =   [(-1,1) (0,1)(-1,0),(0,0),

          (0,1),(1,1),(0,0),(1,0),

          (-1,0),(0,0)(-1,-1),(0,-1)

          (0,0),(1,0)(0,-1),(1,-1)]

我使用以下方法为任意维度 d 生成网格：

grid = [np.linspace(-1.0 , 1.0, num = K+1) for i in range(d)]
res_to_unpack = np.meshgrid(*grid,indexing = 'ij')

哪个有输出：

[array([[-1., -1., -1.],
   [ 0.,  0.,  0.],
   [ 1.,  1.,  1.]]), array([[-1.,  0.,  1.],
   [-1.,  0.,  1.],
   [-1.,  0.,  1.]])]

所以我希望能够为给定的 D 维网格生成上述单元格容器。在给定的 K 上拆分，它是 2 的幂。

我需要这个容器，因此对于每个单元格，我需要引用所有相关的 2^D 点并计算与原点的距离。

编辑澄清

K 应该将网格划分为 K ** D 个单元，其中 (K+1) ** D 个点。每个单元格应有 2 ** D 点数。每个“单元”的体积为 (2/K)^D。

所以对于 K = 4，D = 2

Cells = [ (-1,1),(-0.5,1),(-1,0.5),(-0.5,0.5),
          (-0.5,1),(-0.5,0.5)(0.0,1.0),(0,0.5),
            ...
          (0.0,-0.5),(0.5,-0.5),(0.0,-1.0),(0.5,-1.0),
          (0.5,-1.0),(0.5,-1.0),(1.0,-0.5),(1.0,-1.0)]

这是 TopLeft、TopLeft + Right Over、Bottom Left、Bottom Left + Over Left 的输出。这个集合中有 16 个单元格，每个单元格都有四个坐标。为了增加 K，假设 K = 8。将有 64 个单元格，每个单元格有四个点。

【问题讨论】：

这和your other question有什么不同？

我想我会更具体。

【参考方案1】：

这应该可以满足您的需求：

from itertools import product
import numpy as np

def splitcubes(K, d):
    coords = [np.linspace(-1.0 , 1.0, num=K + 1) for i in range(d)]
    grid = np.stack(np.meshgrid(*coords)).T

    ks = list(range(1, K))
    for slices in product(*([[slice(b,e) for b,e in zip([None] + ks, [k+1 for k in ks] + [None])]]*d)):
        yield grid[slices]

def cubesets(K, d):
    if (K & (K - 1)) or K < 2:
        raise ValueError('K must be a positive power of 2. K: %s' % K)

    return [set(tuple(p.tolist()) for p in c.reshape(-1, d)) for c in splitcubes(K, d)]

2D案例演示

这里有一个 2D 案例的小演示：

import matplotlib.pyplot as plt

def assemblecube(c, spread=.03):
    c = np.array(list(c))
    c = c[np.lexsort(c.T[::-1])]

    d = int(np.log2(c.size))
    for i in range(d):
        c[2**i:2**i + 2] = c[2**i + 1:2**i - 1:-1]

    # get the point farthest from the origin
    sp = c[np.argmax((c**2).sum(axis=1)**.5)]
    # shift all points a small distance towards that farthest point
    c += sp * .1 #np.copysign(np.ones(sp.size)*spread, sp)

    # create several different orderings of the same points so that matplotlib will draw a closed shape
    return [(np.roll(c, i, axis=1) - (np.roll(c, i, axis=1)[0] - c[0])[None,:]).T for i in range(d)]

fig = plt.figure(figsize=(6,6))
ax = fig.gca()

for i,c in enumerate(cubesets(4, 2)):
    for cdata in assemblecube(c):
        p = ax.plot(*cdata, c='C%d' % (i % 9))

ax.set_aspect('equal', 'box')
fig.show()

输出：

出于可视化目的，立方体已稍微分开（因此它们不会重叠并相互覆盖）。

3D案例演示

3D 案例也是如此：

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111, projection='3d')

for i,c in enumerate(cubesets(2,3)):
    for cdata in assemblecube(c, spread=.05):
        ax.plot(*cdata, c=('C%d' % (i % 9)))

plt.gcf().gca().set_aspect('equal', 'box')
plt.show()

输出：

K=4 的演示

这是与上述相同的 2D 和 3D 演示的输出，但带有 K=4：

【讨论】：

如果答案解决了你的问题，你应该mark it as accepted

我想澄清一点，在 K>2 的情况下，假设 K = 4。输出是一个包含 4 个集合的列表，每个集合包含 9 个点。不过，每个单元格仍然应该只包含四个点。

在那种情况下，我真的不明白你想让K 参数做什么。 K=4 与 K=2 的输出应该如何变化？如果您希望我能够向您展示如何正确实现它，您必须在问题的末尾添加一个K=4 的示例，包括所需的输出。

我发布了一个编辑，非常感谢您的帮助，如果问题足够具体，我深表歉意。尽管我是 python 新手，但我正在尝试编辑您的解决方案，所以我很挣扎。

我更新了答案以反映您问题中K 定义的更新。