[toc] 三角剖分有很多种，在刷题中常见的三角剖分一般是 Delaunay三角剖分
这个三角剖分有很多很好的性质

# Delaunay 三角剖分

什么事三角剖分，即对任意多边形进行三角形化且三角形化后的外边缘是这些三角形集合的凸包

# 定义

唯一性
- DT 三角剖分是唯一的，且任意四点不共圆以及外接圆内不会有其他任何一个剖分中的点。
最大化最小角
- 比如上图，在算法中会选择右边作为最终的三角剖分，以保证最小角是最大的 (交换对角线不会导致最小角变大)。

# 性质

最接近：以最接近的三点形成三角形，且个线段均不相交
唯一性：不论从区域何处开始构建，最终都会得到一个结果
最优性：最大化最小角
最规则：如果将剖分结果内的三角形进行排布，最后得到的排列中 DT 剖分的最大
区域性：新增 / 删除 / 移动某一个顶点只会影响周围的点。
具有凸多边形的外壳：三角形外界会形成一个凸多边形

# 算法

OI Wiki 上说的是 O (nlgn) 的算法

算法概念

LL-edge : 左侧子点集的边 RR-edge : 右侧子点集的边 LR-edge : 连接左右剖分产生的新的边 在剖分过程中，可能会删除 LL-edge 或 RR-edge, 在合并时不会增加 LL-edge

先对点集按照 x 坐标进行升序排序
不断分治，不断将点集划分为两个部分，直到其中一个部分点集不超过 3，这个点集就可以立即剖分为一个三角形或线段
左右两个剖分可以合并成一个剖分
合并左右两个剖分：
- 5.1 插入 base LR-edge, 这个是最底部不与任何 LL-edge 以及 RR-edge 相交的 LR-edge
- 5.2 确定下一条 LR-edge, 下一条 LR-edge 可能端点在与 LR-edge 右端点相连的右侧剖分的 RR-edge 上 (6,7,9 点)
- 5.3 找到了新的 LR 边后，判断是否有线相交了，有的话，删除

# 检测可能的端点的标准

对应点对于 RR-edge 和 base_LR-edge 的夹角小于 180 度
base LR-edge 两端点和这个可能点三点构成的圆不包含其他任何可能点

如上图所示，6 号点的外接圆 (绿色) 包含了 9 号点，所以抛弃，选择 7 号点 重复以上步骤，最终就可以得到三角剖分

# 代码

#include <iostream>
#include <algorithm>
#include <cmath>
#include <cstring>
#include <list>
#include <utility>
#include <vector>

const double EPS = 1e-8;
const int MAXV = 10000;

struct Point {
    double x, y;
    int id;

    Point(double a = 0, double b = 0, int c = -1) : x(a), y(b), id(c) {}

    bool operator<(const Point& a) const {
        return x < a.x  (fabs(x - a.x) < EPS && y < a.y);
    }

    bool operator==(const Point& a) const {
        return fabs(x-a.x) < EPS && fabs(y - a.y) < EPS;
    }

    double dist2(const Point& b) {
        return (x - b.x) * (x - b.x) + (y - b.y) * (y - b.y);
    }
};

struct Point3D {
    double x, y, z;

    Point3D(double a = 0, double b = 0, double c = 0) : x(a), y(b), z(c) {}
    //z = 点到原点距离(?)
    Point3D(const Point& p) { x = p.x, y = p.y, z = p.x * p.x + p.y * p.y; }

    Point3D operator-(const Point3D& a) const {
        return Point3D(x - a.x, y - a.y, z - a.z);
    }

    double dot(const Point3D& a) { return x * a.x + y * a.y + z * a.z; }
};

struct Edge {
    int id;
    std::list<Edge>::iterator c;

    Edge(int id = 0) { this->id = id; }
};
//判断v是0还是大于0还是小于0
int cmp(double v) { return fabs(v) > EPS ? (v > 0 ? 1 : -1) : 0; }

//三点叉积
//oa叉积ob
double cross(const Point& o, const Point& a, const Point& b) {
    return (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x);
}

Point3D cross(const Point3D& a, const Point3D& b) {
    return Point3D(a.y * b.z - a.z * b.y, -a.x * b.z + a.z * b.x,
        a.x * b.y - a.y * b.x);
}

int inCircle(const Point& a, Point b, Point c, const Point& p) {
    //先保证叉积向上
    if (cross(a, b, c) < 0) std::swap(b, c);
    Point3D a3(a), b3(b), c3(c), p3(p);
    b3 = b3 - a3, c3 = c3 - a3, p3 = p3 - a3;
    Point3D f = cross(b3, c3);
    //检查是否同向,in: < 0,out > 0
    return cmp(p3.dot(f));
}

//相交检测
//基于如下定理: 
//如果AB叉积AC与AB叉积AD的符号相反,证明CD必然在AB两侧
//同理,如果CD叉积CA与CD叉积CB的符号相反,证明AB必然在CD两侧
//如果都在两侧,则证明这两条线断AB,CD必定相交(不需要跨立)
int intersection(const Point& a, const Point& b, const Point& c,
    const Point& d) {
    return cmp(cross(a, c, b)) * cmp(cross(a, b, d)) > 0 &&
        cmp(cross(c, a, d)) * cmp(cross(c, d, b)) > 0;
}

class Delaunay {
public:
    std::list<Edge> head[MAXV];
    Point p[MAXV];
    int n, rename[MAXV];

    void init(int n, Point p[]) {
        memcpy(this->p, p, sizeof(Point) * n);
        std::sort(this->p, this->p + n);
        for (int i = 0; i < n; i++) rename[p[i].id] = i;
        this->n = n;
        divide(0, n - 1);
    }

    void addEdge(int u, int v) {
        head[u].push_front(Edge(v));
        head[v].push_front(Edge(u));
        head[u].begin()->c = head[v].begin();
        head[v].begin()->c = head[u].begin();
    }

    void divide(int l, int r) {
        //如果点集少于等于三个,则进行第一次三角剖分
        if (r - l <= 2) {
            for (int i = l; i <= r; ++i) {
                for (int j = i + 1; j <= r; ++j) {
                    addEdge(i, j);
                }
            }
            return;
        }

        int mid = (l + r) / 2;
        divide(l, mid);
        divide(mid + 1, r);

        std::list<Edge>::iterator it;
        int nowl = l, nowr = r;

        for (int update = 1; update;) {
            // find left and right convex, lower common tangent
            //找到最低的LR
            update = 0;
            Point ptL = p[nowl], ptR = p[nowr];
            for (it = head[nowl].begin(); it != head[nowl].end(); ++it) {
                Point t = p[it->id];
                //ptr_ptl,ptr_t判定是否t在ptr_ptl下方,是的话,或者共线,则以t为left往下找
                double v = cross(ptR, ptL, t);
                //判定两直线是否相交(也不像跨立实验啊?判断点在直线的某侧?
                if (cmp(v) > 0  (cmp(v) == 0 && ptR.dist2(t) < ptR.dist2(ptL))) {
                    nowl = it->id, update = 1;
                    break;
                }
            }
            if (update) continue;
            //找到了最低的left,找最低的r
            for (it = head[nowr].begin(); it != head[nowr].end(); ++it) {
                Point t = p[it->id];
                double v = cross(ptL, ptR, t);
                if (cmp(v) < 0  (cmp(v) == 0 && ptL.dist2(t) < ptL.dist2(ptR))) {
                    nowr = it->id, update = 1;
                    break;
                }
            }
        }
        //找到了最低的LR 加LR边
        addEdge(nowl, nowr);  // add tangent

        //逐个点判断是否能和LR组成一个圈,且圈内没有其他点
        //需要判定的点是LR边上的几个点
        //找到了新的LR边后,判断是否有线相交了,有的话,删除
        for (int update = 1; true;) {
            update = 0;
            Point ptL = p[nowl], ptR = p[nowr];
            int ch = -1, side = 0;
            for (it = head[nowl].begin(); it != head[nowl].end(); it++) {
                if (cmp(cross(ptL, ptR, p[it->id])) > 0 &&
                    (ch == -1  inCircle(ptL, ptR, p[ch], p[it->id]) < 0)) {
                    ch = it->id, side = -1;
                }
            }
            for (it = head[nowr].begin(); it != head[nowr].end(); it++) {
                if (cmp(cross(ptR, p[it->id], ptL)) > 0 &&
                    (ch == -1  inCircle(ptL, ptR, p[ch], p[it->id]) < 0)) {
                    ch = it->id, side = 1;
                }
            }

            if (ch == -1) break;
            if (side == -1) {
                for (it = head[nowl].begin(); it != head[nowl].end();) {
                    if (intersection(ptL, p[it->id], ptR, p[ch])) {
                        head[it->id].erase(it->c);
                        head[nowl].erase(it++);
                    }
                    else {
                        it++;
                    }
                }
                nowl = ch;
                addEdge(nowl, nowr);
            }
            else {
                for (it = head[nowr].begin(); it != head[nowr].end();) {
                    if (intersection(ptR, p[it->id], ptL, p[ch])) {
                        head[it->id].erase(it->c);
                        head[nowr].erase(it++);
                    }
                    else {
                        it++;
                    }
                }
                nowr = ch;
                addEdge(nowl, nowr);
            }
        }
    }

    std::vector<std::pair<int, int> > getEdge() {
        std::vector<std::pair<int, int> > ret;
        ret.reserve(n);
        std::list<Edge>::iterator it;
        for (int i = 0; i < n; i++) {
            for (it = head[i].begin(); it != head[i].end(); it++) {
                if (it->id < i) continue;
                ret.push_back(std::make_pair(p[i].id, p[it->id].id));
            }
        }
        return ret;
    }

};

//三角剖分

int main()
{
    std::cout << "Hello World!\n";
}

# 后记

由上面我们可以看出，实际上 DT 三角剖分算法是针对有点时的剖分算法，但是在实际情况中我们只有一个完整的多边形，但是在图形渲染流水线中有曲面细分这个流程，这个时候怎么进行三角剖分呢

可以在外部添加一个 AABB 包围盒
对多边形边进行等分 (即辅助点)
以多边形为内部点，AABB 为凸包，进行一次 DT 三角剖分
删除外部的 AABB 和对应的边，得到一个没有内部点的 DT 三角剖分
想办法合理的添加内部点 (如内部三角形的重心等)