计算机视觉：三角剖分算法

算法解析

参考论文Triangulate得到对应的伪代码

subroutine triangulate
input : vertex list
output : triangle list
   initialize the triangle list
   determine the supertriangle
   add supertriangle vertices to the end of the vertex list
   add the supertriangle to the triangle list
   for each sample point in the vertex list
      initialize the edge buffer
      for each triangle currently in the triangle list
         calculate the triangle circumcircle center and radius
         if the point lies in the triangle circumcircle then
            add the three triangle edges to the edge buffer
            remove the triangle from the triangle list
         endif
      endfor
      delete all doubly specified edges from the edge buffer
         this leaves the edges of the enclosing polygon only
      add to the triangle list all triangles formed between the point 
         and the edges of the enclosing polygon
   endfor
   remove any triangles from the triangle list that use the supertriangle vertices
   remove the supertriangle vertices from the vertex list
end

同时，参考了这篇博客对算法进行了优化，提高了运行速度（代码中的method 1标注；method 2为伪代码的实现）

算法实现

// 使用到的自定义的类
// 二维图像上的一个点
class Point{
public:
    double x, y;
    Point(double a, double b) {
        x = a;
        y = b;
    }
    Point(){
        x = 0;
        y = 0;
    }
    bool operator==(const Point& other) {
        return x == other.x && y == other.y;
    }
};
// 二维图像上的一个点对
struct PointPair {
    Point s, t;
    PointPair(Point a, Point b){
        s.x = a.x;
        s.y = a.y;
        t.x = b.x;
        t.y = b.y;
    }
};
// 二维图像上的边，用点对表示
typedef PointPair Edge;
// 二维图像上的一个三角形
struct Triangle {
    Point a,b,c,center;
    double r;
    Triangle(Point p1, Point p2, Point p3, int mode = 0) {
    	// 计算外接圆半径（利用向量法求外接圆圆心和半径）
        if(mode == 1) {
            double A = pow(p1.x, 2) + pow(p1.y, 2);
            double B = pow(p2.x, 2) + pow(p2.y, 2);
            double C = pow(p3.x, 2) + pow(p3.y, 2);
            double G = (p3.y-p2.y) * p1.x + (p1.y-p3.y) * p2.x + (p2.y-p1.y) * p3.x;
            center.x = ((B-C) * p1.y + (C-A) * p2.y + (A-B) * p3.y) / (2*G);
            center.y = ((C-B) * p1.x + (A-C) * p2.x + (B-A) * p3.x) / (2*G);
            r = sqrt(pow(p1.x-center.x, 2) + pow(p1.y-center.y, 2));
        }
        else {
            center.x = 0;
            center.y = 0;
            r = 0;
        }
        a.x = p1.x;
        a.y = p1.y;
        b.x = p2.x;
        b.y = p2.y;
        c.x = p3.x;
        c.y = p3.y;
    }
};

// 三角剖分算法
void Delaunay() {
    vector<Point> vertex;
    vector<Triangle> temp_triangles, triangles;
    double minX, maxX, minY, maxY;
    double maxTriS, maxTriH;
    unsigned char green[3] = {0,255,0};
    unsigned char red[3] = {255,0,0};

    for(int i = 0; i < pointPairList.size(); i++) {
        vertex.push_back(pointPairList[i].s);
    }
	// 初始化一个大三角形，将所有标记点包含起来
    sort(vertex.begin(), vertex.end(), [](const Point& a, const Point& b) { return a.y < b.y; });
    minY = vertex[0].y;
    maxY = vertex[vertex.size()-1].y;
    sort(vertex.begin(), vertex.end(), [](const Point& a, const Point& b) { return a.x < b.x; });
    minX = vertex[0].x;
    maxX = vertex[vertex.size()-1].x;

    maxTriS = (maxX - minX)/2;
    maxTriH = (maxY - minY)*2;
    Point right(maxX+maxTriS+15, maxY+5);
    Point left(minX-maxTriS-15, maxY+5);
    Point top(minX+maxTriS, maxY-maxTriH);
    Triangle maxTriangle(right, left, top, 1);
    temp_triangles.push_back(maxTriangle);
    triangles.push_back(maxTriangle);
	// 三角剖分算法核心代码
    for(int i = 0; i < vertex.size(); i++) {
        vector<Edge> buff, buff_res;
        for(vector<Triangle>::iterator p = temp_triangles.begin(); p != temp_triangles.end(); ){
            // method 1: faster
            if(getDistance(vertex[i], p->center) > p->r && vertex[i].x - p->center.x > p->r) {
                triangles.push_back(*p);
                p = temp_triangles.erase(p);
                continue;
            }
            else if(getDistance(vertex[i], p->center) > p->r && vertex[i].x - p->center.x <= p->r) {
                p++;
                continue;
            }
            else {
                buff.push_back(Edge(p->a, p->b));
                buff.push_back(Edge(p->b, p->c));
                buff.push_back(Edge(p->a, p->c));
                p = temp_triangles.erase(p);
            }
            // method 2: slower
            /*if(getDistance(vertex[i], p->center) < p->r) {
                buff.push_back(Edge(p->a, p->b));
                buff.push_back(Edge(p->b, p->c));
                buff.push_back(Edge(p->a, p->c));
                p = temp_triangles.erase(p);
            } else {
                p++;
            }*/
        }
        int* pos = new int[buff.size()];
        memset(pos, 0, sizeof(int) * buff.size());
        for(int j = 0; j < buff.size()-1; j++) {
            if(pos[j] != 0) continue;
            for(int k = j+1; k < buff.size(); k++) {
                if((buff[j].s == buff[k].s && buff[j].t == buff[k].t) ||
                   (buff[j].s == buff[k].t && buff[j].t == buff[k].s)){
                    pos[k] = 1;
                    pos[j] = 1;
                }
            }
        }
        for(int j = 0; j < buff.size(); j++){
            if(pos[j] == 0) {
                buff_res.push_back(buff[j]);
            }
        }
        for(int j = 0; j < buff_res.size(); j++) {
            Triangle t(buff_res[j].s, buff_res[j].t, vertex[i], 1);
            temp_triangles.push_back(t);
        }
        delete[] pos;
    }
    // 将上述得到的临时三角形存下
    for(int i = 0; i < temp_triangles.size(); i++) {
        triangles.push_back(temp_triangles[i]);
    }
    // 删去所有与外部三角形的端点相关的三角形
    for(vector<Triangle>::iterator p = triangles.begin(); p != triangles.end(); p++) {
        if((*p).a == right || (*p).a == left || (*p).a == top ||
           (*p).b == right || (*p).b == left || (*p).b == top ||
           (*p).c == right || (*p).c == left || (*p).c == top
       ){
           continue;
       }
       else {
           res_triangle.push_back((*p));
       }
    }
}