3D Delaunay三角剖分异常：生成额外面的问题排查与解决

我使用Python 3.11通过scipy.spatial.Delaunay处理点云时，生成的三角剖分结果在曲面四边边界处出现额外面。

点云3D散点图代码 #

import plotly.graph_objects as go
fig = go.Figure()
fig.add_scatter3d(x = puntos[:,0], y = puntos[:,1], z = puntos[:,2],mode='markers', marker=dict(
                size=1,
                color='rgb(0,0,0)',
                opacity=0.8
            ))
fig.update_layout(scene = dict(aspectmode = 'data'))

fig.show()

三角剖分代码 #

import numpy as np
import pandas as pd
from scipy.spatial import Delaunay
import plotly.figure_factory as ff

puntos = pd.read_csv('puntos.csv')
puntos = puntos[['0', '1', '2']]
tri = Delaunay(np.array([puntos[:,0], puntos[:,1]]).T)
simplices = tri.simplices
fig = ff.create_trisurf(x=puntos[:,0], y=puntos[:,1], z=puntos[:,2],
                         simplices=simplices, aspectratio=dict(x=1, y=1, z=0.3))
fig.show()

问题原因 #

1. 2D投影剖分的局限性：你调用Delaunay时仅传入了x、y坐标的2D投影，算法会在平面上对所有点做无约束三角剖分，包括边界外围的点对组合。但原数据是3D曲面，这些平面上的外围三角面在3D空间中会连接曲面边界两端本不相邻的点，形成跨边界的额外面。
2. 无边界约束的剖分逻辑：标准Delaunay三角剖分默认会填充所有点集凸包范围内的区域，不会识别你的点云是有边界的曲面，因此会自动生成边界外的多余三角面。

解决方法 #

方法1：通过法向量过滤额外面 #

计算每个三角面的法向量，与曲面整体法向量对比，剔除方向相反的面：

import numpy as np
import pandas as pd
from scipy.spatial import Delaunay
import plotly.figure_factory as ff

puntos = pd.read_csv('puntos.csv')
puntos = puntos[['0', '1', '2']].values  # 转为numpy数组

# 2D Delaunay剖分
tri = Delaunay(puntos[:, :2])
simplices = tri.simplices

# 过滤方向异常的三角面
valid_simplices = []
for simplex in simplices:
    p1, p2, p3 = puntos[simplex]
    # 计算边向量与法向量
    v1 = p2 - p1
    v2 = p3 - p1
    normal = np.cross(v1, v2)
    # 保留法向量z分量为正的面（根据你的曲面朝向调整）
    if normal[2] > 0:
        valid_simplices.append(simplex)

# 绘制过滤后的三角面
fig = ff.create_trisurf(x=puntos[:,0], y=puntos[:,1], z=puntos[:,2],
                         simplices=valid_simplices, aspectratio=dict(x=1, y=1, z=0.3))
fig.show()

方法2：添加边界约束的2D剖分 #

如果点云边界规则，手动定义边界约束，让算法仅在内部生成三角面：

import numpy as np
import pandas as pd
from scipy.spatial import Delaunay
import plotly.figure_factory as ff

puntos = pd.read_csv('puntos.csv')
puntos = puntos[['0', '1', '2']].values

# 提取边界点（示例：取x、y极值点，需根据实际点云调整）
x_min, x_max = puntos[:,0].min(), puntos[:,0].max()
y_min, y_max = puntos[:,1].min(), puntos[:,1].max()
boundary_mask = ((puntos[:,0] == x_min) | (puntos[:,0] == x_max)) & ((puntos[:,1] == y_min) | (puntos[:,1] == y_max))
boundary_points = puntos[boundary_mask]
# 按顺时针顺序排列边界点
boundary_points = boundary_points[np.lexsort((boundary_points[:,1], boundary_points[:,0]))]

# 生成约束边（相邻点连接，首尾闭合）
constraints = []
for i in range(len(boundary_points)):
    idx1 = np.where((puntos[:,0] == boundary_points[i,0]) & (puntos[:,1] == boundary_points[i,1]))[0][0]
    idx2 = np.where((puntos[:,0] == boundary_points[(i+1)%len(boundary_points),0]) & (puntos[:,1] == boundary_points[(i+1)%len(boundary_points),1]))[0][0]
    constraints.append([idx1, idx2])

# 带约束的Delaunay剖分（需scipy>=1.6.0）
tri = Delaunay(puntos[:, :2], constraints=constraints)
simplices = tri.simplices

fig = ff.create_trisurf(x=puntos[:,0], y=puntos[:,1], z=puntos[:,2],
                         simplices=simplices, aspectratio=dict(x=1, y=1, z=0.3))
fig.show()

方法3：3D Delaunay剖分提取表面 #

直接做3D剖分，提取凸包表面作为有效三角面：

import numpy as np
import pandas as pd
from scipy.spatial import Delaunay
import plotly.figure_factory as ff

puntos = pd.read_csv('puntos.csv')
puntos = puntos[['0', '1', '2']].values

# 3D Delaunay剖分
tri = Delaunay(puntos)
# 提取凸包表面的三角面
surface_faces = tri.convex_hull.tolist()

# 绘制表面
fig = ff.create_trisurf(x=puntos[:,0], y=puntos[:,1], z=puntos[:,2],
                         simplices=surface_faces, aspectratio=dict(x=1, y=1, z=0.3))
fig.show()

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