Python 3中曲面三角化与插值及交线值提取问题

针对你只有5个散点的情况，interp2d确实不适用（它要求输入是规则网格数据），我们可以用径向基函数(RBF)插值来生成曲面，再结合三角剖分和平面相交计算来实现目标。以下是分步实现代码：

1. 依赖库安装 #

先确保安装所需库：

pip install numpy scipy matplotlib

2. 完整代码实现 #

import numpy as np
from scipy.interpolate import Rbf
from scipy.spatial import Delaunay
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 你的原始数据
x = np.array([0.0, 17.67599999997765, 49.08499999996275, 90.57299999985844, 136.60500000044703])
y = np.array([0.0, 45.22349889159747, 66.50303846438841, 114.04427618243405, 187.7707039612985])
z = np.array([0.0, 1.8700000000000045, 1.9539999999999509, 1.3929999999999154, 1.6299999999999955])

# ---------------------- 步骤1：生成插值曲面 ----------------------
# 创建x-y网格（分辨率可调整，这里用20x20）
xi = np.linspace(x.min(), x.max(), 20)
yi = np.linspace(y.min(), y.max(), 20)
xi, yi = np.meshgrid(xi, yi)

# 用RBF插值生成z值（适合散点插值，点少也能稳定生成曲面）
rbf = Rbf(x, y, z, function='multiquadric')
zi = rbf(xi, yi)

# ---------------------- 步骤2：曲面三角化 ----------------------
# 将网格点转为散点格式，用于Delaunay三角剖分
points = np.vstack((xi.flatten(), yi.flatten())).T
tri = Delaunay(points)

# ---------------------- 步骤3：提取与指定平面的交线 ----------------------
# 示例：指定平面为 z = 1.5（你可根据需求修改为其他平面方程，比如ax+by+cz=d）
target_z = 1.5
intersection_points = []

# 遍历每个三角形，检查是否与目标平面相交
for simplex in tri.simplices:
    # 获取三角形三个顶点的z值
    z_vals = zi.flatten()[simplex]
    # 判断三角形是否跨目标平面（有顶点在平面上方，有在下方）
    if (np.any(z_vals > target_z) and np.any(z_vals < target_z)):
        # 遍历三角形的每条边，找与平面的交点
        for i in range(3):
            j = (i + 1) % 3
            z1, z2 = z_vals[i], z_vals[j]
            if (z1 - target_z) * (z2 - target_z) < 0:
                # 线性插值计算交点坐标
                t = (target_z - z1) / (z2 - z1)
                x_inter = points[simplex[i], 0] + t * (points[simplex[j], 0] - points[simplex[i], 0])
                y_inter = points[simplex[i], 1] + t * (points[simplex[j], 1] - points[simplex[i], 1])
                intersection_points.append([x_inter, y_inter, target_z])

# 转成numpy数组方便后续处理
intersection_points = np.array(intersection_points)

# ---------------------- 可视化结果 ----------------------
fig = plt.figure(figsize=(12, 6))

# 子图1：生成的曲面与原始点
ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_trisurf(points[:,0], points[:,1], zi.flatten(), triangles=tri.simplices, alpha=0.5, cmap='viridis')
ax1.scatter(x, y, z, color='red', s=50, label='原始数据点')
ax1.set_xlabel('X')
ax1.set_ylabel('Y')
ax1.set_zlabel('Z')
ax1.legend()
ax1.set_title('生成的曲面')

# 子图2：曲面与指定平面的交线
ax2 = fig.add_subplot(122, projection='3d')
ax2.plot_trisurf(points[:,0], points[:,1], zi.flatten(), triangles=tri.simplices, alpha=0.3, cmap='viridis')
ax2.scatter(intersection_points[:,0], intersection_points[:,1], intersection_points[:,2], color='blue', s=30, label='交线点')
ax2.plot(intersection_points[:,0], intersection_points[:,1], intersection_points[:,2], color='blue', linewidth=2)
ax2.set_xlabel('X')
ax2.set_ylabel('Y')
ax2.set_zlabel('Z')
ax2.legend()
ax2.set_title(f'曲面与z={target_z}平面的交线')

plt.tight_layout()
plt.show()

关键说明 #

• RBF插值：相对于interp2d，它更适合散点数据，即使只有少量点也能生成连续曲面，你可以尝试修改function参数（比如'linear'、'gaussian'）调整曲面平滑度。
• 三角化：用Delaunay对插值后的网格点做三角剖分，得到曲面的三角面片，方便后续计算平面相交。
• 平面相交计算：遍历每个三角面片，判断是否与目标平面交叉，再对交叉的边做线性插值得到交点，最终所有交点组成交线。
• 自定义平面：如果你的目标平面不是z=常数，而是形如a*x + b*y + c*z = d，只需要修改交点判断逻辑，将z_vals替换为a*points[simplex,0] + b*points[simplex,1] + c*z_vals，判断是否跨d即可。

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