《Learning OpenCV 3》Delaunay三角剖分和Voronoi图讲解

OpenCV2、OpenCV3包含三角剖分的接口，但是参考文档里并未介绍，给学习带来了麻烦。

有一本经典的书《学习OpenCV》对其做了详细介绍。但苦于这本书的新版，迟迟没有翻译成中文，所以现在有关OpenCV三角剖分的资料都是关于OpenCV1的，是C语言的接口。

我边学习边分享一下，自己的理解水平有限，如有错漏，谢谢指正。由于时间有限，不逐字翻译了，根据内容来。

原文参考《Learning OpenCV 3》附录A：平面剖分。第923-937页。

---

前面很大一部分是平面剖分的理论，不做翻译，只简单提一下。讲了Delaunay三角剖分的特性和相关的计算方法。Delaunay是一个点集三角化的标准，在所有剖分方案中，满足所有三角形的最小角的和是最大的。Delaunay剖分是唯一的，是最接近规则化的三角网。Delaunay剖分有很多算法实现。OpenCV应用的是逐点插入法，这从函数接口可以得出。

Delaunay三角剖分和Voronoi图是对偶的，这意味着计算出了Delaunay剖分，Voronoi图也就确定了。直观感受一下：

---

重点来了，下面是OpenCV有关函数接口的理解和使用方法说明。

创建Delaunay 或者 Voronoi 剖分

首先需要在内存中开辟一块地方来存储Delaunay剖分。我们也需要一个方框（记住，为了加速计算，算法处理过程中，在这个方框的外面需要一个虚拟的外围三角形）。

为了尽快开始，就假设这些点必须在一个600*600的图像中吧。

// STRUCTURE FOR DELAUNAY SUBDIVISION
//
...
cv::Rect rect(0, 0, 600, 600); // Our outer bounding box
cv::Subdiv2D subdiv(rect); // Create the initial subdivision

这些代码创建了初始的剖分，一个三角形包含一个特定的矩形框。

接下来，我们需要知道怎么插入点。这些点必须是32位float类型的，或者是带有整数坐标值的点（cv::Point）。在后面的案例中，它们会自动转换为float类型。插入点使用cv::Subdiv2D::insert()函数。

（译者注：方便起见，后面很多相关函数省略cv::Subdiv2D命名空间了）

cv::Point2f fp; //This is our point holder
for( int i = 0; i < as_many_points_as_you_want; i++ ) {
// However you want to set points
//
fp = your_32f_point_list[i];
subdiv.insert(fp);
}

现在，点已经输入完毕，我们能够得到Delaunay三角剖分。从Delaunay三角剖分中计算三角形，使用getTiangleList()函数。

vector<cv::Vec6f> triangles;
subdiv.getTriangleList(triangles);

调用之后，，在三角形中的每个Vec6f包含三个顶点：（$x_1,y_1,x_2,y_2,x_3,y_3,$）。

得到对应的Voronoi图使用函数getVoronoiFaceList。

vector<vector<cv::Point2f> > facets;
vector<cv::Point2f> centers;
subdiv.getVoronoiFacetList(vector<int>(), facets, centers);

facets包含Voronoi小面块（译者注：里面的点数据只包括多边形的顶点），centers包含对应的区域中心。

值得一提的是，Delaunay三角化是迭代构建的，意味着，每插入一个点，三角剖分都会更新，所以它是总是更新的。然而，Voronoi图是当你调用calcVoronoi()一次性构建的。可选的是，你可以调用前面提到的getVoronoiFaceList()（它内部调用了calcVoronoi()）来随之更新。

既然，我们已经创建好了一个二维点集的Delaunay剖分以及对应的Voronoi图。下一步就是学习怎么遍历这个剖分。

---

遍历Delaunay剖分

平面剖分基本的数据元素是边，通过序号访问边。通过这个序号还可访问相邻边，附加的参数还可指定想要访问的边同当前边的位置关系。每个边两个端点叫做origin和destination。一个边会和其它边共享这些点。最后，存在一个对应（对偶）边，每个Delaunay剖分的边都有Voronoi剖分的边相对应。

记住一点，在cv::Subdiv2D接口中，对待边总是直接的，这实际上是为了方便。

还有，边有方向。两个点包含两条边，因为区分origin和destination。

根据边访问点

不论Delaunay剖分，还是Voronoi剖分，边都有起点和终点。访问边的端点如下：

int cv::Subdiv2D::edgeOrg( int edge, cv::Point2f* orgpt = 0 ) const;
int cv::Subdiv2D::edgeDst( int edge, cv::Point2f* dstpt = 0 ) const;

edge是输入，是边的序号。参数表第二项返回点本身。函数返回点的序号。

给定点序号，可以得到点的坐标，和相关的边。

cv::Point2f cv::Subdiv2D::getVertex( int vertex, int* firstEdge = 0 ) const;

（译者注：这里的firstEdge和点的关系不清楚，读者可以考证一下）

需要注意，和边一样，点有序号。当然点也有坐标。Subdiv2D接口故意设计成这样，在绝大多数的接口函数中，你主要使用的是边和点的序号。

在剖分中定位点

一个可能发生的情况是：你有一个特定点的位置信息，但是想找到它在剖分中的序号。
相似的情况：可能这个点实际上并不是剖分中的顶点，但是你想找到包含这个点的三角形或小面块。方法locate()把一个点作为输入，返回这个点所在的一条边。或者包含这个点的三角形或面块的一条边（如果这条边不是顶点）。注意，在这种情况下，返回的不一定是距离最近的边，只是简单的返回包含点的三角形或面块的其中一条边。当点是顶点时，locate()也会返回顶点的ID。

int cv::Subdiv2D::locate(
cv::Point2f pt,
int& edge,
int& vertex
);

函数返回值，告诉我们，点的着落位置

cv::Subdiv2D::PTLOC_INSIDE
点落在面块内部，*edge是其中一条边。

cv::Subdiv2D::PTLOC_ON_EDGE
点落在边上 *edge包含这条边。

cv::Subdiv2D::PTLOC_VERTEX
点落在剖分的顶点上， *vertex 包含顶点指针。

cv::Subdiv2D::PTLOC_OUTSIDE_RECT
点落在参考矩形外面，返回指针无效。

cv::Subdiv2D::PTLOC_ERROR
输入参数无效访问

环绕顶点遍历

给定一条边，你可能想访问跟这条边的起点或终点连接的新边。实现这项工作的方法是，我们指定一个起始边，我们绕着它的“头”点或者“尾”点，逆时针，亦或者顺时针搜寻下一条边。这种设计的说明见下图。我们通过函数getEdge()来实现。

int cv:Subdiv2D::getEdge(
int edge,
int nextEdgeType // see text below
) const;

当调用这个函数时，我们需提供当前边和nextEdgeType参数，可选的参数值如下：

• cv::Subdiv2D::NEXT_AROUND_ORG, 绕起点下一条边 (eOnext)
• cv::Subdiv2D::NEXT_AROUND_DST, 绕终点下一条边 (eDnext)
• cv::Subdiv2D::PREV_AROUND_ORG, 绕起点上一条边 (反向 eRnext)
• cv::Subdiv2D::PREV_AROUND_DST, 绕终点上一条边 (反向 eLnext)

怎么遍历完全取决于你，也可以绕三角形或面块遍历，参数值如下：

• cv::Subdiv2D::NEXT_AROUND_LEFT, 环绕左侧面块的下一条边 (eLnext)
• cv::Subdiv2D::NEXT_AROUND_RIGHT,环绕右侧面块的下一条边 (eRnext)
• cv::Subdiv2D::PREV_AROUND_LEFT, 环绕左侧面块的上一条边 (reversed eOnext)
• cv::Subdiv2D::PREV_AROUND_RIGHT, 环绕右侧面块的上一条边(reversed eDnext)

（译者注：“下一条边”的隐含的访问顺序，绕点访问时是逆时针，绕多边形时是逆时针环行）

不用担心是绕Delaunay三角形的边，还是Voronoi图的多边形的边，因为输入edge的序号已经包含这个信息，后面可详细了解边的编号方法。

也可选择方便的调用方式，nextEdge()：

// equivalent to getEdge(edge, cv::Subdiv2D::NEXT_AROUND_ORG)
//
int cv:Subdiv2D::nextEdge(
int edge
) const;

它等价于getEdge()函数按cv::Subdiv2D::NEXT_AROUND_ORG方式调用。当我们想访问环绕一个点的所有边时，这个函数很方便。对一些应用场景很有帮助，比如，从虚拟外接三角形内的某个顶点出发，寻找凸包。

旋转边

假设你手头上有一个边的序号。无论你是从其他函数中得到的，还是想轻率地从某个特定的序号开始遍历整个图，调用下面的函数你可以从Delaunay剖分的边上跳到对应的Voronoi剖分的边上。

int cv::Subdiv2D::rotateEdge(
int edge,
int rotate // get other edges in the same quad-edge: modulo 4 operation
) const;

参数rotate指定了你想旋转的方式，可以选择下列参数指定下一条边，参考下图更易理解：

• 0, 输入的边 (e 下图e 即是)
• 1, 旋转后的边 (eRot)
• 2, 反向边 (反向 e)
• 3, 反向旋转后的边 (反向 eRot)

（译者注：从图中可以看出，旋转边默认是绕起始点逆时针）

关于顶点和边更多的知识

顶点及其命名顺序

由于Delaunay剖分初始化的方式，下面的事实总是成立的：

1. 0号顶点是空顶点，没有坐标。
2. 接下来，1、2、3号顶点是给定外围矩形外侧的“虚拟”顶点，每个都被设置成距离点集很远的位置。
   （译者注：理解为一个虚无缥缈的位置，不确定的位置，我将这3个点看作外接三角形的3个顶点，因为三角形也是虚构的。）
3. 接下来的点都是点集上的点，提供给Subdiv2D对象。

边及其命名顺序

在Subdiv2D对象里的每条边都被赋予一个整数值，这些整数被4个一组使用，每4个号码代表的边是相关联的：

• edge % 4 == 0
  一条 Delaunay 边

• edge % 4 == 1
  垂直于初始边的Voronoi 边

• edge % 4 == 2
  和初始边方向相反

• edge % 4 == 3
  上面Voronoi 边的反向

虚拟边和空边

0号边是空边，不指向任何地方（或者，更准确地说，它的两个端点是0号顶点-也是空的）。

1、2、3号边总是连接虚拟顶点的虚拟Delaunay边。指未固定的虚拟边，因为它的两个顶点都是虚拟的。（译者注：因为边的两头都是虚无缥缈，边当然也是了，而且没有一端是实打实的点。）

空边的起点和终点都是（0，0）。

旋转空边的结果，会得到另一个空边。从空点开始的“第一条边”也是一个空边，随后用nextEdge() 产生的边也一样。

从任何虚拟顶点访问的“第一条边”总是连接到另一个虚拟顶点。（译者注：“第一条边”怎么理解？）

确定外接三角形

既然，当我们对一个点集创建Delaunay剖分时，前3个点总是构建出一个外接三角形（不包括0号点）。我们可以通过下面的方式访问这三个顶点：

Point2f outer_vtx[3];
for( int i = 0; i < 3; i++ ) {
outer_vtx[i] = subdiv.getVertex(i+1);
}

我们也能得到外接三角形的3条边：

int outer_edges[3];
outer_edges[0] = 1*4;
outer_edges[1] = subdiv.getEdge(outer_edges[0], Subdiv2D::NEXT_AROUND_LEFT);
outer_edges[2] = subdiv.getEdge(outer_edges[1], Subdiv2D::NEXT_AROUND_LEFT);

（译者注：这么说来，4号边总是外接三角形的一条边）

确定凸包，在凸包上环绕访问

回忆一下，根据构造函数Subdiv2D(rect)，我们用一个外接矩形初始化了Delaunay剖分。基于此，下面的叙述成立：

• 如果有这样一条边，它的起点和终点都在矩形外侧，然后这条边在剖分的虚构外接三角形内。这样的边，我们叫做未固定的虚拟边。

• 如果有这样一条边，它的两个端点分布在矩形内外两侧，然后内侧的点在点集的凸包上。凸包上的每个点都和虚构的外围三角形的两个顶点相连，而且这两条边的序号是相邻的。我们把这样一端连在矩形内，一端连在矩形外虚拟点上的边，叫做固定的虚拟边。

基于以上事实，我们可以快速找到凸包。例如：从顶点1、2、3开始，我们知道这3个点是虚构外接三角形上的3个虚拟顶点。我们可以使用nextEdge() 可以迅速产生所有的固定的虚拟边的集合（简单的拒绝未固定的虚拟边）。然后调用rotateEdge()，取反向边，然后再调用1次，或者2次nextEdge,就会落在凸包的边上。准确的讲，一条固定的虚拟边对应一条凸包的边，这些边的集合就是凸包。

（译者注：
解决一些疑惑
1. nextEdge()可访问点四周的所有边，所以从1，2，3开始会得到全部的固定的虚拟边。
2. 通过边可以获得起点和终点的编号，通过边两端点的序号可以区分未固定的虚拟边、固定的虚拟边、未固定的虚拟边。
3. rotateEdge(2)将边的起点从虚构三角形顶点上转移到凸包的顶点上。
4. 调用1次或2次nextEdge()恰恰验证了凸包的顶点跟虚拟三角形有两条连线。
）

使用示例

我们可以使用locate()环绕Delaunay三角形的边逐步访问。在下面的例子中，写了一个函数实现给定一个点，在包含点的Delaunay三角形的每条边上做些什么事：

void locate_point(
cv::Subdiv2D& subdiv,
const cv::Point2f& fp,
...
) {
    int e;
    int e0 = 0;
    int vertex = 0;
    subdiv.locate( fp, e0, vertex );
    if( e0 > 0 ) {
        e = e0;
        do // Always 3 edges -- this is a triangulation, after all.
        {
        // [Insert your code here]
        //
        // Do something with e ...
        e = subdiv.getEdge( e, cv::Subdiv2D::NEXT_AROUND_LEFT );
        }
        while( e != e0 );
    }
}

给定一个点，我们也可以调用如下函数找到最近的点：

int Subdiv2D::findNearest(
cv::Point2f pt,
cv::Point2f* nearestPt
);

跟locate()不同，findNearest()会返回剖分中距离最近的顶点的整数ID。输入的点不必落在面块或三角形内。值得注意的是，这个函数不是const函数，因为当没有更新数据时它会计算Voronoi图。

类似的，我们可以环绕Voronoi面块访问，然后画出来它。

void draw_subdiv_facet(
cv::Mat& img,
cv::Subdiv2D& subdiv,
int edge
) {
    int t = edge;
    int i, count = 0;
    vector<cv::Point> buf;
    // Count number of edges in facet
    do{
        count++;
        t = subdiv.getEdge( t, cv::Subdiv2D::NEXT_AROUND_LEFT );
    } while (t != edge );
    // Gather points
    //
    buf.resize(count);
    t = edge;
    for( i = 0; i < count; i++ ) {
        cv::Point2f pt;
        if( subdiv.edgeOrg(t, &pt) <= 0 )
            break;
        buf[i] = cv::Point(cvRound(pt.x), cvRound(pt.y));
        t = subdiv.getEdge( t, cv::Subdiv2D::NEXT_AROUND_LEFT );
    }
    // Around we go
    //
    if( i == count ){
        cv::Point2f pt;
        subdiv.edgeDst(subdiv.rotateEdge(edge, 1), &pt);
        fillConvexPoly(
        img, buf,
        cv::Scalar(rand()&255,rand()&255,rand()&255),
        8, 0
        );
    vector< vector<cv::Point> > outline;
    outline.push_back(buf);
    polylines(img, outline, true, cv::Scalar(), 1, cv::LINE_AA, 0);
    draw_subdiv_point( img, pt, cv::Scalar(0,0,0) );
    }
}