Skip to main content
 首页 » 编程设计

python之将 Python MeshGrid 拆分为单元格

2025年05月04日20kevingrace

问题陈述

需要将 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 上拆分,该 K 是 2 的幂。

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

编辑澄清

K 应该将网格划分为 K ** D 个单元格,具有 (K+1) ** D 个点。每个 Cell 应该有 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 个单元格,每个单元格有四个点。

请您参考如下方法:

这应该给你你所需要的:

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)] 

二维案例演示

这是二维案例的一个小演示:

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: