如何从单元中心坐标计算八叉树计算网格的尺寸

问题背景 #

现有八叉树计算网格：蓝色网格为边长1的立方体，红色网格边长0.5，绿色网格边长0.25。仅通过单元中心坐标cellcenters.csv自动计算每个网格的尺寸并可视化，但原有Python代码输出结果错误。

原错误代码 #

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import cKDTree

file_path = 'cellcenters.csv'
cell_centers_df = pd.read_csv(file_path, header=None, names=['x', 'y', 'z'])

def plot_cell_boundaries(ax, boundaries, color='b', label=None):
    for boundary in boundaries:
        x_min, x_max, y_min, y_max, z_min, z_max = boundary
    
        ax.plot([x_min, x_max], [y_min, y_min], [z_min, z_min], color=color)
        ax.plot([x_min, x_max], [y_max, y_max], [z_min, z_min], color=color)
        ax.plot([x_min, x_max], [y_min, y_min], [z_max, z_max], color=color)
        ax.plot([x_min, x_max], [y_max, y_max], [z_max, z_max], color=color)

        ax.plot([x_min, x_min], [y_min, y_max], [z_min, z_min], color=color)
        ax.plot([x_max, x_max], [y_min, y_max], [z_min, z_min], color=color)
        ax.plot([x_min, x_min], [y_min, y_max], [z_max, z_max], color=color)
        ax.plot([x_max, x_max], [y_min, y_max], [z_max, z_max], color=color)

        ax.plot([x_min, x_min], [y_min, y_min], [z_min, z_max], color=color)
        ax.plot([x_max, x_max], [y_min, y_min], [z_min, z_max], color=color)
        ax.plot([x_min, x_min], [y_max, y_max], [z_min, z_max], color=color)
        ax.plot([x_max, x_max], [y_max, y_max], [z_min, z_max], color=color)

kdtree = cKDTree(cell_centers_df[['x', 'y', 'z']])
distances, indices = kdtree.query(cell_centers_df[['x', 'y', 'z']], k=2)

cell_sizes = distances[:, 1] / 2

def calculate_variable_size_cell_boundaries(df, sizes):
    cell_boundaries = []

    for i, row in df.iterrows():
        x_center, y_center, z_center = row['x'], row['y'], row['z']
        cell_size = sizes[i]
    
        x_min, x_max = x_center - cell_size, x_center + cell_size
        y_min, y_max = y_center - cell_size, y_center + cell_size
        z_min, z_max = z_center - cell_size, z_center + cell_size
    
        cell_boundaries.append([x_min, x_max, y_min, y_max, z_min, z_max])

    return np.array(cell_boundaries)

variable_cell_boundaries = calculate_variable_size_cell_boundaries(cell_centers_df, cell_sizes)

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

plot_cell_boundaries(ax, variable_cell_boundaries, color='b')

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('Cell Boundaries (Wrong)')

plt.show()

错误原因解析 #

原代码通过最近邻距离的一半计算网格半长，这是错误的：

• 八叉树网格是层级化的，不同层级的网格尺寸差异大，大网格的最近邻可能是小网格，导致计算出的尺寸远小于实际值。
• 例如，大网格中心(1.5,0.5,1.5)的最近邻是小网格中心(1.25,0.25,1.25)，计算出的半长约为0.216，远小于实际的0.5（对应边长1的网格）。

修正方案：利用八叉树网格的层级特性 #

八叉树网格的每个单元格尺寸是基础尺寸的2^k倍（这里基础尺寸为0.25），且单元格中心的小数部分始终是对应尺寸半长的奇数倍：

• 边长1的网格：半长0.5，中心小数部分为0.5（如0.5,1.5）
• 边长0.5的网格：半长0.25，中心小数部分为0.25或0.75（如1.25,1.75）
• 边长0.25的网格：半长0.125，中心小数部分为0.125/0.375/0.625/0.875（如1.125,1.375）

修正后代码 #

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

file_path = 'cellcenters.csv'
cell_centers_df = pd.read_csv(file_path, header=None, names=['x', 'y', 'z'])

def plot_cell_boundaries(ax, boundaries, color='b', label=None):
    for boundary in boundaries:
        x_min, x_max, y_min, y_max, z_min, z_max = boundary
    
        # 绘制立方体的12条边
        ax.plot([x_min, x_max], [y_min, y_min], [z_min, z_min], color=color)
        ax.plot([x_min, x_max], [y_max, y_max], [z_min, z_min], color=color)
        ax.plot([x_min, x_max], [y_min, y_min], [z_max, z_max], color=color)
        ax.plot([x_min, x_max], [y_max, y_max], [z_max, z_max], color=color)

        ax.plot([x_min, x_min], [y_min, y_max], [z_min, z_min], color=color)
        ax.plot([x_max, x_max], [y_min, y_max], [z_min, z_min], color=color)
        ax.plot([x_min, x_min], [y_min, y_max], [z_max, z_max], color=color)
        ax.plot([x_max, x_max], [y_min, y_max], [z_max, z_max], color=color)

        ax.plot([x_min, x_min], [y_min, y_min], [z_min, z_max], color=color)
        ax.plot([x_max, x_max], [y_min, y_min], [z_min, z_max], color=color)
        ax.plot([x_min, x_min], [y_max, y_max], [z_min, z_max], color=color)
        ax.plot([x_max, x_max], [y_max, y_max], [z_min, z_max], color=color)

def get_cell_half_size(coord):
    """根据坐标的小数部分判断网格半长"""
    frac = coord - np.floor(coord)
    # 检查属于哪个层级
    if np.isclose(frac, 0.5):
        return 0.5
    elif np.isclose(frac, 0.25) or np.isclose(frac, 0.75):
        return 0.25
    elif any(np.isclose(frac, val) for val in [0.125, 0.375, 0.625, 0.875]):
        return 0.125
    else:
        # 处理异常情况，默认最小尺寸
        return 0.125

def calculate_octree_cell_boundaries(df):
    cell_boundaries = []
    cell_colors = []
    for _, row in df.iterrows():
        x, y, z = row['x'], row['y'], row['z']
        # 取三个维度中最大的半长（八叉树网格三个维度尺寸一致）
        half_sizes = [get_cell_half_size(x), get_cell_half_size(y), get_cell_half_size(z)]
        half_size = max(half_sizes)
        
        x_min, x_max = x - half_size, x + half_size
        y_min, y_max = y - half_size, y + half_size
        z_min, z_max = z - half_size, z + half_size
        
        cell_boundaries.append([x_min, x_max, y_min, y_max, z_min, z_max])
        # 按尺寸分配颜色
        if half_size == 0.5:
            cell_colors.append('blue')
        elif half_size == 0.25:
            cell_colors.append('red')
        else:
            cell_colors.append('green')
    return np.array(cell_boundaries), cell_colors

# 计算边界和颜色
cell_boundaries, cell_colors = calculate_octree_cell_boundaries(cell_centers_df)

# 可视化
fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111, projection='3d')

# 按颜色绘制每个网格
for boundary, color in zip(cell_boundaries, cell_colors):
    plot_cell_boundaries(ax, [boundary], color=color)

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('Octree Cell Boundaries (Corrected)')
ax.legend(['Size 1 (Blue)', 'Size 0.5 (Red)', 'Size 0.25 (Green)'])

plt.show()

代码说明 #

1. get_cell_half_size函数：通过坐标的小数部分判断该点所属网格的半长，利用八叉树网格的坐标规律实现精准匹配。
2. calculate_octree_cell_boundaries函数：为每个中心计算网格边界，并按尺寸分配对应颜色（蓝/红/绿）。
3. 可视化优化：按网格层级分颜色绘制，更直观区分不同尺寸的八叉树网格。

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