Python中用OpenCV与igraph实现Delaunay三角剖分的异常问题

背景与初始问题 #

使用Python 3.9.13实现实时点集Delaunay三角剖分并分析图属性，流程为：

• 用OpenCV（opencv-python 4.6.0.66）的Subdiv2D生成剖分
• 转换为igraph（igraph 0.10.3）可分析的图

初始问题：每隔几帧igraph生成的图会错乱，仅80%的时间能正常工作，排查发现trilist生成的顶点列表长度存在波动。

初始完整代码 #

import time
import numpy
import cv2
import igraph as ig

# Draw a point
def draw_point(img, p, color):
    cv2.circle(img, (int(p[0]), int(p[1])), 2, color, 0)

# Get a triangulation
def get_delaunay(subdiv):
    return subdiv.getTriangleList()

# Draw delaunay triangles
def draw_delaunay(img, subdiv, delaunay_color):
    triangleList = subdiv.getTriangleList()
    size = img.shape
    r = (0, 0, size[1], size[0])

    for t in triangleList:
        pt1 = (int(t[0]), int(t[1]))
        pt2 = (int(t[2]), int(t[3]))
        pt3 = (int(t[4]), int(t[5]))
        cv2.line(img, pt1, pt2, delaunay_color, 1, cv2.LINE_AA, 0)
        cv2.line(img, pt2, pt3, delaunay_color, 1, cv2.LINE_AA, 0)
        cv2.line(img, pt3, pt1, delaunay_color, 1, cv2.LINE_AA, 0)

if __name__ == '__main__':
    NUM_PART = 500
    SIZE = 1000
    REPEAT = 10
    for iteration in range(REPEAT):
        positions = numpy.random.randint(0, SIZE, size=(NUM_PART, 2))
        print("There is {p} positions. And {up} unique position".format(p=len(positions), up=len(numpy.unique(positions, axis=1))))
        rect = (0, 0, SIZE, SIZE)
        subdiv = cv2.Subdiv2D(rect)
        timer = time.time()
        print("There is {} points in subdiv".format(len(positions)))
        for p in positions:
            p = p.astype("float32")
            subdiv.insert(p)
        trilist = get_delaunay(subdiv)
        print("Took {}s".format(round(time.time() - timer, 12)))
        print("there is {} triangles in trilist".format(len(trilist)))
        opencv_image = numpy.zeros((SIZE, SIZE, 3))
        draw_delaunay(opencv_image, subdiv, (255, 255, 255))
        for p in positions:
            draw_point(opencv_image, p, (0, 0, 255))
        timer = time.time()
        n_vertices = NUM_PART
        g = ig.Graph(n=n_vertices, )
        g.vs["name"] = range(NUM_PART)
        print("graph name vector of length {l}:\n{v}".format(l=len(g.vs["name"]), v=g.vs["name"]))
        positionsx = [SIZE - pos for pos in positions[:, 0]]
        g.vs["x"] = positions[:, 0]
        print("graph x vector of length {l}:\n{v}".format(l=len(g.vs["x"]), v=g.vs["x"]))
        g.vs["y"] = positions[:, 1]
        print("graph y vector of length {l}:\n{v}".format(l=len(g.vs["y"]), v=g.vs["y"]))
        print("Graph took {}s".format(round(time.time() - timer, 12)))
        list_vtx = []
        for tri in trilist:
            vertex1, _ = subdiv.findNearest((tri[0], tri[1]))
            vertex2, _ = subdiv.findNearest((tri[2], tri[3]))
            vertex3, _ = subdiv.findNearest((tri[4], tri[5]))
            list_vtx.extend([vertex3, vertex2, vertex1])
        list_cleared = list(dict.fromkeys(list_vtx))
        list_cleared.sort()
        print("list cleared of length {len}: {lst}".format(len=len(list_cleared), lst=list_cleared))
        for tri in trilist:
            vertex1, _ = subdiv.findNearest((tri[0], tri[1]))
            vertex2, _ = subdiv.findNearest((tri[2], tri[3]))
            vertex3, _ = subdiv.findNearest((tri[4], tri[5]))
            g.add_edges([(vertex1 - 4, vertex2 - 4), (vertex2 - 4, vertex3 - 4), (vertex3 - 4, vertex1 - 4)])
        g.simplify()
        nodes = g.vs.indices
        print(nodes)
        print(subdiv)
        igraph_image = numpy.zeros((SIZE, SIZE, 3))
        for point in g.vs:
            draw_point(igraph_image, (point["x"], point["y"]), (0, 0, 255))
        for edge in g.es:
            cv2.line(igraph_image, (int(g.vs["x"][edge.tuple[0]]), int(g.vs["y"][edge.tuple[0]])), (int(g.vs["x"][edge.tuple[1]]), int(g.vs["y"][edge.tuple[1]])), (255, 255, 255), 1, cv2.LINE_AA, 0)
        numpy_horizontal = numpy.hstack((opencv_image, igraph_image))
        cv2.imshow('L: opencv || R: igraph', numpy_horizontal)
        cv2.waitKey(0)

修改后的代码与新问题 #

按照建议修改代码后，图错乱问题消失，但出现部分边丢失的情况（即使在之前正常的场景下也存在），修改后的核心代码片段如下：

x_pos = [0] * NUM_PART   # create 0-filled array of x-positions
y_pos = [0] * NUM_PART   # create 0-filled array of y-positions
edges = []               # create empty array of edges
for tri in trilist:
    vertex1 = subdiv.findNearest((tri[0], tri[1]))[0] - 4
    vertex2 = subdiv.findNearest((tri[2], tri[3]))[0] - 4
    vertex3 = subdiv.findNearest((tri[4], tri[5]))[0] - 4
    x_pos[vertex1] = tri[0]
    y_pos[vertex1] = tri[1]
    x_pos[vertex2] = tri[2]
    y_pos[vertex2] = tri[3]
    x_pos[vertex3] = tri[4]
    y_pos[vertex3] = tri[5]
    edges.append((vertex1, vertex2))
    edges.append((vertex2, vertex3))
    edges.append((vertex2, vertex3))  # 此处重复添加了同一条边
g = ig.Graph(NUM_PART, edges)
g.vs["name"] = range(NUM_PART)
g.vs["x"] = x_pos
g.vs["y"] = y_pos
g.simplify()

问题原因与解决方法 #

问题原因 #

1. 边添加逻辑错误：三角形的三条边应为(vertex1, vertex2)、(vertex2, vertex3)、(vertex3, vertex1)，但代码中重复添加了(vertex2, vertex3)，遗漏了(vertex3, vertex1)，导致每个三角形少一条边。
2. 重复点未处理：输入的positions可能存在重复点，subdiv.insert会自动去重，此时NUM_PART大于实际顶点数，导致索引不匹配。
3. 硬编码索引偏移风险：依赖-4跳过OpenCV自动添加的4个边界顶点，若顶点存储逻辑变化会直接导致索引错误。

解决方法 #

1. 修正边添加逻辑：将重复的边替换为正确的第三条边：

   edges.append((vertex1, vertex2))
   edges.append((vertex2, vertex3))
   edges.append((vertex3, vertex1))  # 修正为正确的第三条边
   

2. 提前去重顶点：生成positions后先去重，保证插入到Subdiv2D的顶点数与igraph的顶点数一致：

   # 生成positions后执行去重
   positions = numpy.unique(positions, axis=0)
   NUM_PART = len(positions)  # 更新顶点数
   

3. 替换索引获取方式：通过Subdiv2D.getPoints()直接获取顶点列表，跳过前4个边界点，避免findNearest的索引偏差：

   # 获取所有顶点，跳过前4个边界点
   custom_points = subdiv.getPoints()[4:]
   # 建立坐标到索引的映射
   point_to_idx = {tuple(p): i for i, p in enumerate(custom_points)}
   # 遍历三角形时直接通过坐标找索引
   for tri in trilist:
       p1 = (tri[0], tri[1])
       p2 = (tri[2], tri[3])
       p3 = (tri[4], tri[5])
       idx1 = point_to_idx[p1]
       idx2 = point_to_idx[p2]
       idx3 = point_to_idx[p3]
       edges.extend([(idx1, idx2), (idx2, idx3), (idx3, idx1)])
   

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