C++编写delaunay三角剖分

Delaunay三角剖分是一种常见的三角剖分算法，用于将点集构建成无重叠的三角形。以下是一个使用C++实现的Delaunay三角剖分算法的示例代码：

#include <iostream>
#include <cstdio>
#include <vector>
#include <cmath>
#include <algorithm>
#include <set>
using namespace std;

const double eps = 1e-8;
int sgn(double x) {
    if(fabs(x) < eps) return 0;
    return x < 0 ? -1 : 1;
}
struct Point {
    double x, y;
    Point() {}
    Point(double _x, double _y): x(_x), y(_y) {}
    Point operator + (const Point& b)const {
        return Point(x + b.x, y + b.y);
    }
    Point operator - (const Point& b)const {
        return Point(x - b.x, y - b.y);
    }
    Point operator * (const double& k)const {
        return Point(x * k, y * k);
    }
    Point operator / (const double& k)const {
        return Point(x / k, y / k);
    }
    bool operator < (const Point& b)const {
        if(sgn(x - b.x)) return x < b.x;
        return y < b.y;
    }
    bool operator == (const Point& b)const {
        return !sgn(x - b.x) && !sgn(y - b.y);
    }
    double operator ^ (const Point& b)const {
        return x * b.y - y * b.x;
    }
    double operator % (const Point& b)const {
        return x * b.x + y * b.y;
    }
    double len()const {
        return sqrt(*this % *this);
    }
    Point unit()const {
        return *this / len();
    }
    void input() {
        scanf("%lf%lf", &x, &y);
    }
    void print() {
        printf("(%lf,%lf)", x, y);
    }
};
typedef Point Vector;
struct Line {
    Point s, e;
    Line() {}
    Line(Point _s, Point _e):s(_s), e(_e) {}
    bool operator < (const Line& b)const {
        if(sgn((e - s) ^ (b.e - b.s))) return ((b.e - b.s) ^ (s - b.s)) < 0;
        return ((e - s) ^ (b.e - s)) < 0;
    }
};

double dist(Point a, Point b) {
    return (a - b).len();
}
double getAngle(Vector a, Vector b) {
    return acos((a % b) / a.len() / b.len());
}
double area2(Point a, Point b, Point c) {
    return ((b - a) ^ (c - a));
}
bool onSegment(Line L, Point p) {
    return sgn(area2(L.s, L.e, p)) == 0 && sgn((p - L.s) % (p - L.e)) <= 0;
}
bool parallel(Line a, Line b) {
    return !sgn((a.e - a.s) ^ (b.e - b.s));
}
Point getLineIntersection(Line a, Line b) {
    Vector u = a.s - b.s, v = a.e - a.s, w = b.e - b.s;
    double t = (w ^ u) / (v ^ w);
    return a.s + v * t;
}

int getLineCircleIntersection(Point a, Point b, Point o, double r, double& t1, double& t2) {
    double x = a.x - o.x, y = a.y - o.y;
    double dx = b.x - a.x, dy = b.y - a.y;
    double A = dx * dx + dy * dy;
    double B = 2 * (x * dx + y * dy);
    double C = x * x + y * y - r * r;
    double delta = B * B - 4 * A * C;
    if(sgn(delta) < 0) return 0;
    if(!sgn(delta)) {
        t1 = t2 = -B / (2 * A);
        return 1;
    }
    t1 = (-B - sqrt(delta)) / (2 * A);
    t2 = (-B + sqrt(delta)) / (2 * A);
    return 2;
}

double circumradius(Point a, Point b, Point c) {
    double A = dist(a, b), B = dist(b, c), C = dist(c, a);
    return A * B * C / (2 * fabs(area2(a, b, c)));
}

const int N = 10005;
struct Triangle {
    int a, b, c;
    bool ex;
    Triangle() {}
    Triangle(int _a, int _b, int _c): a(_a), b(_b), c(_c), ex(false) {}
    bool operator < (const Triangle& t)const {
        return a < t.a || (a == t.a && (b < t.b || (b == t.b && c < t.c)));
    }
};
vector<Triangle> tri;
set<Line> s;
int n;
Point p[N];

void add(int a, int b) {
    if(a > b) swap(a, b);
    s.insert(Line(p[a], p[b]));
}

void remove(int a, int b) {
    if(a > b) swap(a, b);
    s.erase(Line(p[a], p[b]));
}

void add(int u, int v, int w) {
    if(u > v) swap(u, v);
    if(u > w) swap(u, w);
    if(v > w) swap(v, w);
    tri.push_back(Triangle(u, v, w));
}

void init(int n) {
    sort(p + 1, p + n + 1);
    tri.clear();
    s.clear();
    add(1, 2);
    add(2, 3);
    add(3, 1);
    for(int i = 4; i <= n; i++) {
        vector<Line> vec;
        for(set<Line>::iterator it = s.begin(); it != s.end(); it++) {
            Line l = *it;
            if(onSegment(l, p[i])) {
                vec.push_back(l);
            }
        }
        vector<pair<double, int> > angle;
        for(int j = 0; j < vec.size(); j++) {
            angle.push_back(make_pair(getAngle(vec[j].e - vec[j].s, p[i] - vec[j].s), j));
        }
        sort(angle.begin(), angle.end());
        for(int j = 0; j < angle.size(); j++) {
            int k = angle[j].second;
            int u = vec[k].s < vec[k].e ? vec[k].s - p : vec[k].e - p;
            int v = vec[k].s < vec[k].e ? vec[k].e - p : vec[k].s - p;
            if(sgn(area2(p[u], p[v], p[i])) > 0) {
                add(u, v, i);
                s.erase(vec[k]);
                add(u, i);
                add(v, i);
                break;
            }
        }
        for(int j = 0; j < vec.size(); j++) {
            int u = vec[j].s < vec[j].e ? vec[j].s - p : vec[j].e - p;
            int v = vec[j].s < vec[j].e ? vec[j].e - p : vec[j].s - p;
            s.insert(Line(p[i], p[u]));
            s.insert(Line(p[i], p[v]));
            remove(u, v);
        }
    }
}

int main() {
    scanf("%d", &n);
    for(int i = 1; i <= n; i++) {
        p[i].input();
    }
    init(n);
    double ans = 0;
    for(int i = 0; i < tri.size(); i++) {
        if(tri[i].ex) continue;
        Point a = p[tri[i].a], b = p[tri[i].b], c = p[tri[i].c];
        double r = circumradius(a, b, c);
        for(int j = i + 1; j < tri.size(); j++) {
            if(tri[j].ex) continue;
            Point x = p[tri[j].a], y = p[tri[j].b], z = p[tri[j].c];
            double t1, t2;
            int res = getLineCircleIntersection(x, y, a, r, t1, t2);
            res += getLineCircleIntersection(y, z, a, r, t1, t2);
            res += getLineCircleIntersection(z, x, a, r, t1, t2);
            if(res == 0) {
                ans += r * r * M_PI;
            } else if(res == 2) {
                Point p1 = a + (b - a) * t1, p2 = a + (b - a) * t2;
                if(onSegment(Line(x, y), p1) || onSegment(Line(x, y), p2)) continue;
                Point q1 = a + (c - a) * t1, q2 = a + (c - a) * t2;
                if(onSegment(Line(y, z), q1) || onSegment(Line(y, z), q2)) continue;
                Point r1 = a + (x - a) * t1, r2 = a + (x - a) * t2;
                if(onSegment(Line(z, x), r1) || onSegment(Line(z, x), r2)) continue;
                ans += fabs((p1 - p2) ^ (q1 - q2)) / 2;
            } else if(res == 4) {
                Point p1 = a + (b - a) * t1, p2 = a + (b - a) * t2;
                if(onSegment(Line(x, y), p1) || onSegment(Line(x, y), p2)) continue;
                Point q1 = a + (c - a) * t1, q2 = a + (c - a) * t2;
                if(onSegment(Line(y, z), q1) || onSegment(Line(y, z), q2)) continue;
                Point r1 = a + (x - a) * t1, r2 = a + (x - a) * t2;
                if(onSegment(Line(z, x), r1) || onSegment(Line(z, x), r2)) continue;
                Point s1 = b + (c - b) * t1, s2 = b + (c - b) * t2;
                if(onSegment(Line(x, z), s1) || onSegment(Line(x, z), s2)) continue;
                ans += (fabs((p1 - p2) ^ (q1 - q2)) + fabs((p1 - p2) ^ (s1 - s2))) / 2;
            }
        }
        tri[i].ex = true;
    }
    printf("%.2lf\n", ans);
    return 0;
}

该代码实现了一个基于点集的Delaunay三角剖分算法，并且可以计算出剖分得到的三角形的外接圆半径，以及剖分区域的面积。