如何求CloughTocher2DInterpolator插值曲面与平面的交线？

问题背景 #

给定一组稀疏点集，通过CloughTocher2DInterpolator生成连续曲面，需求解该曲面与自定义平面z = a*x + b*y + c的交线。现有尝试中，三角剖分求交效果差，curve_fit因插值器在凸包外返回NaN导致优化卡顿。

原始插值代码：

import numpy as np
from scipy.interpolate import CloughTocher2DInterpolator

rng = np.random.default_rng()
x = rng.random(10) - 0.5
y = rng.random(10) - 0.5
z = np.hypot(x, y)
X = np.linspace(min(x), max(x))
Y = np.linspace(min(y), max(y))
X, Y = np.meshgrid(X, Y)  # 2D grid for interpolation
interp = CloughTocher2DInterpolator(list(zip(x, y)), z)

可行解决方案 #

方案1：网格等值线提取（基于CloughTocher插值器） #

利用密集网格采样插值曲面，通过等值线提取工具直接获取交线，实现简单且效果稳定。

步骤：

1. 生成覆盖原始点凸包的密集网格，计算插值曲面的Z值；
2. 计算插值曲面与目标平面的差值，屏蔽凸包外的NaN区域；
3. 提取差值为0的等值线，即为交线。

代码示例：

import numpy as np
from scipy.interpolate import CloughTocher2DInterpolator
import matplotlib.pyplot as plt

# 生成原始点集与插值器（同问题背景代码）
rng = np.random.default_rng()
x = rng.random(10) - 0.5
y = rng.random(10) - 0.5
z = np.hypot(x, y)
interp = CloughTocher2DInterpolator(list(zip(x, y)), z)

# 自定义平面参数
a, b, c = 1, 1, 0.1

# 生成密集网格（密度越高精度越高）
X = np.linspace(min(x), max(x), 200)
Y = np.linspace(min(y), max(y), 200)
X, Y = np.meshgrid(X, Y)
Z_interp = interp(X, Y)

# 计算曲面与平面的差值
diff = Z_interp - (a * X + b * Y + c)
diff[np.isnan(diff)] = np.nan  # 屏蔽凸包外无效区域

# 提取等值线（值为0的线段即为交线）
contours = plt.contour(X, Y, diff, levels=[0])
plt.close()  # 仅提取数据，不显示绘图窗口

# 收集交线点集
intersection_points = []
for contour in contours.collections:
    for path in contour.get_paths():
        intersection_points.extend(path.vertices)

intersection_points = np.array(intersection_points)
print("交线点集：\n", intersection_points)

优缺点：实现简单，依赖成熟工具；精度由网格密度决定，密度越高计算量越大。

方案2：RBF插值器+根查找（支持凸包外交线） #

更换为Rbf径向基函数插值器，其可在凸包外正常插值（无NaN），结合根查找方法求解交线。

步骤：

1. 用Rbf构建插值器（可选cubic核函数保证连续光滑）；
2. 定义差值函数f(x,y) = 插值曲面z - 平面z；
3. 遍历x范围，用根查找求解对应y值，收集交线点。

代码示例：

import numpy as np
from scipy.interpolate import Rbf
from scipy.optimize import root

# 生成原始点集
rng = np.random.default_rng()
x = rng.random(10) - 0.5
y = rng.random(10) - 0.5
z = np.hypot(x, y)

# 创建RBF插值器（cubic核保证C2连续）
interp_rbf = Rbf(x, y, z, function='cubic')

# 自定义平面参数
a, b, c = 1, 1, 0.1

# 定义待求解的差值函数
def equation(vars):
    x_val, y_val = vars
    return interp_rbf(x_val, y_val) - (a * x_val + b * y_val + c)

# 遍历x范围求解对应y值
x_range = np.linspace(min(x)-0.1, max(x)+0.1, 50)
intersection_points = []
for x_val in x_range:
    # 初始猜测值基于平面公式推导
    y_guess = (interp_rbf(x_val, 0) - a*x_val - c) / b
    try:
        sol = root(equation, [x_val, y_guess])
        if sol.success:
            intersection_points.append(sol.x)
    except Exception:
        continue

intersection_points = np.array(intersection_points)
print("交线点集：\n", intersection_points)

优缺点：支持凸包外交线求解；但计算量大于CloughTocher，核函数选择会影响插值效果。

方案3：三角剖分逐面片求交（高精度高效） #

利用CloughTocher2DInterpolator的Delaunay三角剖分结构，仅在可能与平面相交的三角形内采样求交，精度与效率平衡。

步骤：

1. 获取插值器的三角剖分结构；
2. 遍历每个三角形，判断是否与平面相交（顶点差值符号变化）；
3. 在相交三角形内采样，通过线性插值找到精确零点，收集交线段。

代码示例：

import numpy as np
from scipy.interpolate import CloughTocher2DInterpolator
from scipy.spatial import Delaunay

# 生成原始点集与插值器
rng = np.random.default_rng()
x = rng.random(10) - 0.5
y = rng.random(10) - 0.5
z = np.hypot(x, y)
points = np.column_stack((x, y))
tri = Delaunay(points)
interp = CloughTocher2DInterpolator(points, z)

# 自定义平面参数
a, b, c = 1, 1, 0.1
plane_z = lambda x_val, y_val: a*x_val + b*y_val + c

intersection_segments = []

# 遍历每个三角形面片
for simplex in tri.simplices:
    # 获取三角形顶点与对应z值
    p1, p2, p3 = points[simplex]
    z1, z2, z3 = z[simplex]
    
    # 计算顶点处曲面与平面的差值
    d1 = z1 - plane_z(p1[0], p1[1])
    d2 = z2 - plane_z(p2[0], p2[1])
    d3 = z3 - plane_z(p3[0], p3[1])
    
    # 判断是否相交：差值有正负变化或存在零点
    if not ((d1*d2 < 0) or (d2*d3 <0) or (d3*d1 <0) or any(np.isclose([d1,d2,d3], 0))):
        continue
    
    # 在三角形内生成采样点
    u = np.linspace(0, 1, 20)
    v = np.linspace(0, 1, 20)
    u, v = np.meshgrid(u, v)
    mask = u + v <= 1
    u, v = u[mask], v[mask]
    # 转换为三角形内的x,y坐标
    x_tri = p1[0]*(1-u-v) + p2[0]*u + p3[0]*v
    y_tri = p1[1]*(1-u-v) + p2[1]*u + p3[1]*v
    z_tri = interp(x_tri, y_tri)
    diff_tri = z_tri - plane_z(x_tri, y_tri)
    
    # 查找符号变化的点，插值得到精确零点
    sign_changes = np.where(np.diff(np.sign(diff_tri)) != 0)[0]
    for idx in sign_changes:
        # 线性插值求零点
        t = -diff_tri[idx] / (diff_tri[idx+1] - diff_tri[idx])
        x_intersect = x_tri[idx] + t*(x_tri[idx+1]-x_tri[idx])
        y_intersect = y_tri[idx] + t*(y_tri[idx+1]-y_tri[idx])
        intersection_segments.append([x_intersect, y_intersect])

# 去重处理
intersection_segments = np.unique(np.round(intersection_segments, decimals=6), axis=0)
print("交线点集：\n", intersection_segments)

优缺点：仅在相交区域计算，效率高且精度可控；实现稍复杂，需处理三角形内采样与零点插值。

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