带周期性边界条件的3D数组分块及重构问题求助

分割函数（支持任意N和n） #

核心思路是利用数组周期性扩展处理边界跨越问题，按n//2的步长生成子立方体，确保相邻子立方体重叠一半体积，同时覆盖整个原数组（包括边界跨越的情况）。

import numpy as np

def split_into_overlapping_subcubes(arr, n):
    N = arr.shape[0]
    stride = n // 2
    # 生成2倍大小的周期性扩展数组，方便直接切片获取跨边界的子立方体
    periodic_arr = np.tile(arr, (2, 2, 2))
    
    # 生成每个维度的起始索引，确保覆盖整个原数组
    dim_starts = np.arange(0, N, stride)
    # 补充最后一个起始点，保证末尾区域被覆盖（当N无法被stride整除时）
    if dim_starts[-1] + n < N + stride:
        dim_starts = np.append(dim_starts, N - stride)
    
    # 生成所有三维起始坐标的组合
    starts = np.array(np.meshgrid(dim_starts, dim_starts, dim_starts)).T.reshape(-1, 3)
    
    # 提取每个子立方体并添加通道维度
    subcubes = []
    for x, y, z in starts:
        subcube = periodic_arr[x:x+n, y:y+n, z:z+n]
        subcubes.append(subcube[..., np.newaxis])
    
    return np.array(subcubes)

函数说明 #

• 输入：形状为(N,N,N)的3D numpy数组arr，子立方体边长n
• 输出：形状为(N_subcubes, n, n, n, 1)的子立方体数组，其中N_subcubes是子立方体总数（每个维度的起始点数量的三次方）
• 边界处理：通过np.tile将原数组扩展为2倍大小的周期性数组，跨边界的子立方体直接从扩展数组中切片，自动包含原数组另一侧的元素

重构函数（从子立方体还原原数组） #

针对重叠区域的多份数据，采用平均加权的方式还原原数组，同样通过取模处理周期性坐标映射。

def reconstruct_from_subcubes(subcubes, N):
    n = subcubes.shape[1]
    stride = n // 2
    # 初始化结果数组和覆盖计数数组
    result = np.zeros((N, N, N), dtype=subcubes.dtype)
    count = np.zeros((N, N, N), dtype=int)
    
    # 生成与分割时一致的起始坐标
    dim_starts = np.arange(0, N, stride)
    if dim_starts[-1] + n < N + stride:
        dim_starts = np.append(dim_starts, N - stride)
    starts = np.array(np.meshgrid(dim_starts, dim_starts, dim_starts)).T.reshape(-1, 3)
    
    # 遍历每个子立方体，累加对应原数组位置的数值
    for idx, (x_start, y_start, z_start) in enumerate(starts):
        # 计算子立方体每个位置对应的原数组坐标（取模实现周期性映射）
        x_coords = (x_start + np.arange(n)) % N
        y_coords = (y_start + np.arange(n)) % N
        z_coords = (z_start + np.arange(n)) % N
        xx, yy, zz = np.meshgrid(x_coords, y_coords, z_coords, indexing='ij')
        
        # 累加数值并计数覆盖次数
        result[xx, yy, zz] += subcubes[idx, :, :, :, 0]
        count[xx, yy, zz] += 1
    
    # 取平均得到最终重构结果
    result /= count
    return result

函数说明 #

• 输入：形状为(N_subcubes, n, n, n, 1)的子立方体数组subcubes，原数组边长N
• 输出：形状为(N,N,N)的重构数组
• 重叠处理：每个原数组元素的最终值是所有覆盖它的子立方体对应位置的平均值，保证重叠区域的还原精度

测试示例 #

# 生成测试数据
N = 624
n = 128
arr = np.random.rand(N, N, N)

# 分割
subcubes = split_into_overlapping_subcubes(arr, n)
print(f"分割后子立方体数组形状: {subcubes.shape}")

# 重构
reconstructed_arr = reconstruct_from_subcubes(subcubes, N)

# 验证重构精度（误差应接近0，受浮点数精度影响）
print(f"最大重构误差: {np.max(np.abs(arr - reconstructed_arr)):.6f}")

内容的提问来源于stack exchange，提问作者Sam K