基于PyVista利用网格面法线高程生成体素建筑高度图遇问题

• 斜侧面误判：你之前用face_normals[:,2] > 0.0会把所有朝上的面都算进去，包括建筑侧面的斜向面（比如屋檐、斜坡），这些面不属于顶面但法线z分量为正，导致额外计数。
• 同位置多单元：建筑有凸起、阳台这类结构时，同一个(x,y)网格位置会对应多个顶面单元（比如上层阳台和下层屋顶），自然会超过32x32=1024的数量。
• OBJ体素冗余：输入的OBJ文件可能包含重复体素面（比如地板和建筑底部的重叠面），也会增加无效单元。

1. 严格筛选纯顶面单元 #

把法线判断阈值提高，只保留几乎完全朝上的面（z轴法线接近1），排除斜侧面：

upward_face_indices = np.where(face_normals[:, 2] > 0.99)[0]

2. 网格对齐+去重取最高高度 #

因为体素是规则32x32网格，先把单元中心映射到对应网格位置，再对每个位置保留最高的高程值，最终得到标准32x32高度图：

import pyvista as pv
import numpy as np

mesh = pv.read("voxeled_withfloor.obj")

# 提取顶面单元
elv = mesh.elevation(preference="cell")
elv.compute_normals(cell_normals=True)
face_normals = elv.cell_normals
# 只保留近乎垂直朝上的顶面
upward_face_indices = np.where(face_normals[:, 2] > 0.99)[0]
upper_face_mesh = elv.extract_cells(upward_face_indices)

# 获取单元中心的三维坐标
cell_centers = upper_face_mesh.cell_centers()
centers = cell_centers.points  # 形状为(N,3)，对应x,y,z

# 计算32x32网格的步长和边界
x_min, x_max = centers[:,0].min(), centers[:,0].max()
y_min, y_max = centers[:,1].min(), centers[:,1].max()
x_step = (x_max - x_min) / 31  # 32个点对应31个间隔
y_step = (y_max - y_min) / 31

# 将每个中心坐标映射到网格索引
x_indices = np.round((centers[:,0] - x_min) / x_step).astype(int)
y_indices = np.round((centers[:,1] - y_min) / y_step).astype(int)

# 过滤掉超出0-31范围的异常索引
valid_mask = (x_indices >=0) & (x_indices <32) & (y_indices >=0) & (y_indices <32)
x_indices = x_indices[valid_mask]
y_indices = y_indices[valid_mask]
z_values = centers[:,2][valid_mask]

# 生成32x32高度图，每个网格取最高z值
height_map = np.zeros((32,32))
for x, y, z in zip(x_indices, y_indices, z_values):
    if z > height_map[y, x]:  # 索引顺序根据你的坐标系调整，可能需要交换x/y
        height_map[y, x] = z

# 最终得到标准32x32高度图
print(height_map.shape)  # 输出(32,32)

3. 调试验证 #

可以用PyVista绘制高度图的热力图，确认是否匹配建筑顶面形态：

grid = pv.UniformGrid()
grid.dimensions = (33, 33, 2)
grid.spacing = (x_step, y_step, 1)
grid.origin = (x_min, y_min, 0)
grid.cell_data["Height"] = height_map.flatten()

pv.plot(grid, cmap="terrain")

• 网格索引的顺序（x对应列还是行）需要根据你的体素坐标系调整，可能需要交换height_map[y, x]中的x和y。
• 如果还是有冗余单元，先对原始mesh合并重复面：mesh = mesh.merge()

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