循环处理多几何图形时triangle库触发段错误问题求助

我需要将Shapely多边形分解为三角形，约束条件是三角形顶点只能位于多边形边上。

目前代码多数情况运行正常，但处理某个特定多边形时触发段错误：

（进程退出码139，被信号11:SIGSEGV中断）：POLYGON ((545155.9709965222 5378933.039350484,545155.8786686319 5378932.570194424, 545153.9828366871 5378934.823291995, 545154.3704756413 5378934.7733652275, 545155.9709965222 5378933.039350484))

更奇怪的是：逐行执行代码处理该特定多边形时完全正常，但放入循环遍历整个GeoDataFrame时就崩溃。

运行环境：Python 3.12，triangle库版本20230923。

另外求Python替代方案——Shapely自带的三角剖分会生成内部顶点，不满足我的约束需求。

我的代码如下：

import os
import geopandas as gpd
from shapely.geometry import Point, MultiPolygon, Polygon, GeometryCollection
from shapely.validation import explain_validity
import shapely.geometry as geom
import numpy as np
import triangle as tr


def is_triangle(geometry):
    if isinstance(geometry, Polygon):
        # 单个多边形：判断是否有3个顶点（闭合多边形首尾重复，所以坐标数为4）
        return len(list(geometry.exterior.coords)) == 4
    elif isinstance(geometry, MultiPolygon):
        # 多面多边形：检查每个子多边形都是三角形
        return all(len(list(poly.exterior.coords)) == 4 for poly in geometry.geoms)
    else:
        return False  # 不支持的几何类型


def process_geometry(geometry):
    # 检查几何有效性，无效则尝试用buffer(0)修复
    if not geometry.is_valid:
        print(f"Invalid geometry: {explain_validity(geometry)}")
        geometry = geometry.buffer(0)

    if is_triangle(geometry):
        return [geometry]  # 返回列表保持格式一致
    elif isinstance(geometry, Polygon):
        return triangulate_polygon(geometry)  # 剖分并返回三角形列表
    elif isinstance(geometry, MultiPolygon):
        triangles = []
        for poly in geometry.geoms:
            triangles.extend(triangulate_polygon(poly))  # 剖分每个子多边形
        return triangles
    else:
        raise ValueError("Geometry type not supported")

def triangulate_polygon(geometry):
    # 移除近重复顶点，避免浮点精度问题
    cleaned_coords = remove_near_duplicate_vertices(np.array(geometry.exterior.coords))
    polygon = Polygon(cleaned_coords)

    coords = np.array(polygon.exterior.coords)

    # 检查退化或共线段
    for i in range(len(coords) - 1):
        dist = np.linalg.norm(coords[i] - coords[i + 1])
        if dist < 1e-8:
            print(f"Nearly degenerate segment between points {i} and {i + 1}")

    # 确保有足够的唯一点用于剖分
    if len(coords) < 3:
        return []
    segments = [[i, (i + 1) % len(coords)] for i in range(len(coords) - 1)]

    # 检查线段索引是否越界
    for segment in segments:
        if any(idx >= len(coords) for idx in segment):
            print(f"Invalid segment index {segment}")

    # 使用triangle库执行剖分
    triangulation = tr.triangulate({'vertices': coords, 'segments': segments}, 'p')

    triangles = create_triangles_from_triangulation(triangulation)

    return triangles


def create_triangles_from_triangulation(triangulation):
    """将剖分结果转为Shapely多边形对象"""
    triangles = []
    if 'triangles' in triangulation:
        coords = triangulation['vertices']
        for tri_indices in triangulation['triangles']:
            triangle_coords = coords[tri_indices]
            triangles.append(geom.Polygon(triangle_coords))
    return triangles


all_triangles = []

for id, tri in non_triangles_gdf.iterrows():
    print(tri)
    print(id)
    print(tri['geometry'])
    geom = tri['geometry']
    all_triangles.extend(process_geometry(geom))

1. 补全缺失的函数 #

代码中调用了remove_near_duplicate_vertices但未实现，这可能导致循环中出现隐藏的状态异常，先补全该函数：

def remove_near_duplicate_vertices(coords, threshold=1e-8):
    """移除距离小于阈值的相邻顶点，确保多边形闭合"""
    cleaned = [coords[0]]
    for point in coords[1:]:
        if np.linalg.norm(point - cleaned[-1]) > threshold:
            cleaned.append(point)
    # 保证首尾点一致，维持多边形闭合
    if np.linalg.norm(cleaned[0] - cleaned[-1]) > threshold:
        cleaned.append(cleaned[0])
    return np.array(cleaned)

2. 解决triangle库的段错误 #

可能原因： #

• triangle库的C底层在循环中存在内存泄漏或状态残留，导致累积错误
• 特定多边形的顶点过于接近，触发底层的退化处理逻辑崩溃

修复方案： #

• 进程隔离：在循环中用子进程处理每个多边形，避免状态残留。示例：

  from multiprocessing import Pool
  
  def process_single_geom(geom):
      try:
          return process_geometry(geom)
      except Exception as e:
          print(f"Error processing geometry: {e}")
          return []
  
  # 批量处理
  with Pool() as pool:
      results = pool.map(process_single_geom, non_triangles_gdf['geometry'].tolist())
  for res in results:
      all_triangles.extend(res)
  

• 坐标缩放：将特定多边形的坐标放大（比如乘以1e6），剖分后再缩小，避免浮点精度问题：

  def triangulate_polygon(geometry):
      cleaned_coords = remove_near_duplicate_vertices(np.array(geometry.exterior.coords))
      # 缩放坐标
      scale = 1e6
      scaled_coords = cleaned_coords * scale
      polygon = Polygon(scaled_coords)
  
      coords = np.array(polygon.exterior.coords)
      # ... 原有检查逻辑 ...
  
      triangulation = tr.triangulate({'vertices': coords, 'segments': segments}, 'p')
  
      # 缩放回原坐标
      triangles = []
      if 'triangles' in triangulation:
          coords = triangulation['vertices'] / scale
          for tri_indices in triangulation['triangles']:
              triangle_coords = coords[tri_indices]
              triangles.append(geom.Polygon(triangle_coords))
      return triangles
  

1. 使用trimesh库（易安装，支持约束剖分） #

trimesh的平面网格生成会严格保留原多边形的顶点，满足约束需求：

import trimesh
import shapely.geometry as geom

def triangulate_with_trimesh(polygon):
    # 将Shapely多边形转为trimesh路径
    path = trimesh.load_path(polygon.exterior.coords)
    # 生成带约束的平面三角网格
    mesh = path.to_planar_mesh()
    # 转为Shapely三角形列表
    triangles = []
    for face in mesh.faces:
        tri_coords = mesh.vertices[face]
        # 闭合多边形
        tri_coords = np.vstack([tri_coords, tri_coords[0]])
        triangles.append(geom.Polygon(tri_coords))
    return triangles

2. 使用CGAL的约束Delaunay剖分（精度更高） #

通过pycgal调用CGAL的约束三角剖分功能，需要先安装CGAL（可通过conda安装pycgal）：

from CGAL.CGAL_Triangulation_2 import Constrained_Delaunay_triangulation_2
from CGAL.CGAL_Kernel import Point_2
import shapely.geometry as geom

def triangulate_with_cgal(polygon):
    # 提取多边形顶点（去掉闭合的重复点）
    coords = np.array(polygon.exterior.coords)[:-1]
    # 转为CGAL点对象
    cgal_points = [Point_2(x, y) for x, y in coords]
    # 初始化约束Delaunay剖分
    cdt = Constrained_Delaunay_triangulation_2()
    # 添加多边形边作为约束
    for i in range(len(cgal_points)):
        next_idx = (i + 1) % len(cgal_points)
        cdt.insert_constraint(cgal_points[i], cgal_points[next_idx])
    # 提取所有有限面（三角形）
    triangles = []
    for face in cdt.finite_faces():
        # 获取三个顶点坐标
        pts = [(face.vertex(j).point().x(), face.vertex(j).point().y()) for j in range(3)]
        # 闭合多边形
        pts.append(pts[0])
        triangles.append(geom.Polygon(pts))
    return triangles

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