基于VTK/PyVista的等值线点连接问题求助

我正在从零实现等值线算法（不使用任何contouring或interpolation库函数），用PyVista创建自定义网格并渲染。我的iso_contour函数生成了(x, y, 0)格式的等值线点，但无法将这些点正确连接成闭合的等值线区域。

我的代码 #

import pyvista as pv 
import numpy as np

def f(x, y):
    return np.sin(10*x)+np.cos(4*y)-np.cos(3*x*y)

x = np.arange(0, 1, 0.05)
y = np.arange(0, 1, 0.05)
XX, YY = np.meshgrid(x, y)
data = f(XX, YY)

isovalue = 0.15

def iso_contour(array, target_value):
    contours = []
    rows, cols = array.shape
    
    # Loop through each cell in the array
    for i in range(rows-1):
        for j in range(cols-1):
            cell_verts = []
            
            # Check each corner of the cell
            for corner in [(i,j), (i+1,j), (i+1,j+1), (i,j+1)]:
                x, y = corner
                cell_verts.append((x, y, array[x,y]))


                # Check each edge of the cell
            for k in range(4):
                p1, p2 = cell_verts[k], cell_verts[(k+1)%4]
                if (p1[2] >= target_value) != (p2[2] >= target_value):
                    # Calculate intersection point
                    t = (target_value - p1[2]) / (p2[2] - p1[2])
                    x, y = p1[:2] + t * (np.array(p2[:2]) - np.array(p1[:2]))
                    contours.append((y, x, 0.0))
    
    return contours

contour_points = iso_contour(data, isovalue)
pl = pv.Plotter()
pl.add_mesh(grid, show_edges = True)
for i in range(len(contour_points)-1):
    line = pv.Line(contour_points[i], contour_points[i+1])
    pl.add_mesh(line, color='r', line_width=5)
pl.show()

我期望得到闭合连贯的等值线图形，但实际输出是杂乱无章的线段，没有形成正确的闭合区域。请问是否有VTK过滤器或其他方法可以帮我实现预期的等值线形状？

问题根源 #

你的代码只是遍历所有网格单元格，把每个单元格上的等值线交点按顺序收集到一个列表里，但这些点并没有按照等值线的连通性分组——相邻的点可能属于完全不同的等值线段，直接两两连线自然会混乱。

方法一：用VTK过滤器快速整理交点 #

VTK的vtkStripper可以自动识别连通的点并连成折线，步骤如下：

1. 将交点数据转为PyVista的PolyData对象：

points_np = np.array(contour_points)
poly = pv.PolyData(points_np)

2. 使用vtkStripper处理：

from vtkmodules.vtkFiltersCore import vtkStripper

stripper = vtkStripper()
stripper.SetInputData(poly)
stripper.Update()
stripped_poly = pv.wrap(stripper.GetOutput())

3. 渲染修正后的等值线：

# 先补全你代码中未定义的grid
grid = pv.StructuredGrid(XX, YY, np.zeros_like(XX))
grid["data"] = data

pl = pv.Plotter()
pl.add_mesh(grid, show_edges=True)
pl.add_mesh(stripped_poly, color='r', line_width=5)
pl.show()

这个方法能快速修复连线问题，适合不想修改核心算法的场景。

方法二：改进等值线算法，正确分组连通点 #

如果要保持从零实现的思路，需要给每个交点记录边信息，通过匹配邻边来构建连通的等值线：

1. 修改iso_contour函数，返回分组后的等值线点列表：

def iso_contour(array, target_value):
    # 存储交点：(坐标, 边标识)，边标识格式为(i,j,edge_idx)，edge_idx对应单元格的四条边
    contour_info = []
    rows, cols = array.shape
    
    for i in range(rows-1):
        for j in range(cols-1):
            cell_verts = []
            for corner in [(i,j), (i+1,j), (i+1,j+1), (i,j+1)]:
                x, y = corner
                cell_verts.append((x, y, array[x,y]))
            
            for k in range(4):
                p1, p2 = cell_verts[k], cell_verts[(k+1)%4]
                if (p1[2] >= target_value) != (p2[2] >= target_value):
                    t = (target_value - p1[2]) / (p2[2] - p1[2])
                    x, y = p1[:2] + t * (np.array(p2[:2]) - np.array(p1[:2]))
                    # 记录当前交点所属的单元格边
                    edge_id = (i, j, k)
                    contour_info.append( ((y, x, 0.0), edge_id) )
    
    # 按连通性分组交点
    contours = []
    used_indices = set()
    
    for idx, (point, edge_id) in enumerate(contour_info):
        if idx in used_indices:
            continue
        # 开始构建一条等值线
        current_line = [point]
        used_indices.add(idx)
        current_edge = edge_id
        
        while True:
            # 找到当前边的相邻边（比如单元格的下边对应上方单元格的上边）
            i, j, k = current_edge
            if k == 0:
                neighbor_edge = (i-1, j, 2)
            elif k == 1:
                neighbor_edge = (i, j+1, 3)
            elif k == 2:
                neighbor_edge = (i+1, j, 0)
            elif k == 3:
                neighbor_edge = (i, j-1, 1)
            else:
                break
            
            # 查找邻边对应的交点
            found = False
            for next_idx, (next_point, next_edge) in enumerate(contour_info):
                if next_idx not in used_indices and next_edge == neighbor_edge:
                    current_line.append(next_point)
                    used_indices.add(next_idx)
                    current_edge = next_edge
                    found = True
                    break
            if not found:
                # 检查是否闭合，若闭合则添加起点完成闭环
                if np.linalg.norm(np.array(current_line[0]) - np.array(current_line[-1])) < 1e-6:
                    current_line.append(current_line[0])
                break
        contours.append(current_line)
    
    return contours

2. 渲染分组后的等值线：

contour_lines = iso_contour(data, isovalue)
grid = pv.StructuredGrid(XX, YY, np.zeros_like(XX))
grid["data"] = data

pl = pv.Plotter()
pl.add_mesh(grid, show_edges=True)
for line_points in contour_lines:
    line_poly = pv.PolyData(line_points)
    pl.add_mesh(line_poly, color='r', line_width=5)
pl.show()

这个方法从根源上解决了点不连通的问题，生成的等值线完全符合预期。

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