Voronoi 分割算法c++代码

以下是Voronoi分割算法的C ++代码示例：

#include <iostream>
#include <vector>
#include <queue>
#include <cmath>
using namespace std;

struct Point {
    double x, y;
    Point(double _x, double _y) : x(_x), y(_y) {}
};

struct Edge {
    Point start, end;
    Edge(Point _s, Point _e) : start(_s), end(_e) {}
};

struct Event {
    double x;
    Point p;
    int type;
    Event(double _x, Point _p, int _t) : x(_x), p(_p), type(_t) {}
    bool operator<(const Event& other) const {
        return x < other.x;
    }
};

struct Site {
    Point p;
    int sitenum;
    Site(double x, double y, int n) : p(x, y), sitenum(n) {}
};

struct SiteCompare {
    bool operator()(const Site& l, const Site& r) const {
        if (l.p.y < r.p.y) return true;
        else if (l.p.y == r.p.y && l.p.x < r.p.x) return true;
        else return false;
    }
};

struct GraphEdge {
    double x1, y1, x2, y2;
    int site1, site2;
    GraphEdge(double _x1, double _y1, double _x2, double _y2, int s1, int s2) : x1(_x1), y1(_y1), x2(_x2), y2(_y2), site1(s1), site2(s2) {}
};

class Voronoi {
private:
    vector<double> site_x, site_y;
    vector<Site> sites;
    priority_queue<Event> eventqueue;
    vector<Edge> edges;
    vector<Point> vertices;
    vector<int> edgecnt;
    double xmin, xmax, ymin, ymax, deltax, deltay;
    int siteidx;
    vector<vector<int>> neighbors;
    vector<bool> processed;
    vector<GraphEdge> graphedges;

    void push_site(double x, double y) {
        site_x.push_back(x);
        site_y.push_back(y);
        siteidx++;
    }

    void read_sites() {
        double x, y;
        while (cin >> x >> y) {
            push_site(x, y);
        }
    }

    void compute_bounding_box() {
        xmin = site_x[0];
        xmax = site_x[0];
        for (int i = 1; i < siteidx; i++) {
            if (site_x[i] < xmin) xmin = site_x[i];
            if (site_x[i] > xmax) xmax = site_x[i];
        }
        deltax = xmax - xmin;

        ymin = site_y[0];
        ymax = site_y[0];
        for (int i = 1; i < siteidx; i++) {
            if (site_y[i] < ymin) ymin = site_y[i];
            if (site_y[i] > ymax) ymax = site_y[i];
        }
        deltay = ymax - ymin;
    }

    void initialize() {
        siteidx = 0;
        read_sites();
        compute_bounding_box();
        int n = site_x.size();
        for (int i = 0; i < n; i++) {
            sites.push_back(Site(site_x[i], site_y[i], i));
        }
        sort(sites.begin(), sites.end(), SiteCompare());
        processed.resize(n, false);
        neighbors.resize(n);
    }

    void insert_event(double x, double y) {
        Point p(x, y);
        eventqueue.emplace(x, p, 0);
    }

    bool circle(Point a, Point b, Point c, double& x, double& y, double& r) {
        double x1 = a.x;
        double y1 = a.y;
        double x2 = b.x;
        double y2 = b.y;
        double x3 = c.x;
        double y3 = c.y;

        double fabsy1y2 = abs(y1 - y2);
        double fabsy2y3 = abs(y2 - y3);

        if (fabsy1y2 < 1e-6 && fabsy2y3 < 1e-6) return false;

        double xc, yc;
        if (fabsy1y2 < 1e-6) {
            double m2 = -((x3 - x2) / (y3 - y2));
            double mx2 = (x2 + x3) / 2.0;
            double my2 = (y2 + y3) / 2.0;
            xc = (x2 + x1) / 2.0;
            yc = m2 * (xc - mx2) + my2;
        }
        else if (fabsy2y3 < 1e-6) {
            double m1 = -((x2 - x1) / (y2 - y1));
            double mx1 = (x1 + x2) / 2.0;
            double my1 = (y1 + y2) / 2.0;
            xc = (x3 + x2) / 2.0;
            yc = m1 * (xc - mx1) + my1;
        }
        else {
            double m1 = -((x2 - x1) / (y2 - y1));
            double m2 = -((x3 - x2) / (y3 - y2));
            double mx1 = (x1 + x2) / 2.0;
            double mx2 = (x2 + x3) / 2.0;
            double my1 = (y1 + y2) / 2.0;
            double my2 = (y2 + y3) / 2.0;
            xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
            yc = m1 * (xc - mx1) + my1;
        }

        x = xc;
        y = yc;
        r = sqrt(pow(x1 - xc, 2) + pow(y1 - yc, 2));
        return true;
    }

    void process_event() {
        while (!eventqueue.empty()) {
            Event e = eventqueue.top();
            if (e.type == 0) {
                insert_site_event(e);
            }
            else {
                insert_circle_event(e);
            }
            eventqueue.pop();
        }
    }

    void insert_site_event(Event e) {
        Site newsite(e.p.x, e.p.y, siteidx++);
        edges.push_back(Edge(newsite.p, newsite.p));
        int ns = sites.size();
        for (int i = 0; i < ns; i++) {
            double dx = newsite.p.x - sites[i].p.x;
            double dy = newsite.p.y - sites[i].p.y;
            double dist = sqrt(dx*dx + dy*dy);
            edges.push_back(Edge(newsite.p, sites[i].p));
            edges.push_back(Edge(sites[i].p, newsite.p));
            if (dist > 0) {
                if (circle(sites[i].p, newsite.p, Point(e.p.x, e.p.y), edges.back().start.x, edges.back().start.y, edges.back().end.x)) {
                    edges.pop_back();
                    edges.pop_back();
                }
            }
        }
        sites.push_back(newsite);
    }

    void insert_circle_event(Event e) {
        int i = e.arc;
        if (i > 0 && i < siteidx - 1) {
            int left = -1, right = -1;
            if (check_circle_event(i, left, right, e.p.x)) {
                Point lp = sites[left].p;
                Point rp = sites[right].p;
                double x, y, r;
                circle(lp, sites[i].p, rp, x, y, r);
                vertices.push_back(Point(x, y));
                eventqueue.emplace(x, Point(x, y), 1);
                edges.push_back(Edge(lp, rp));
                edgecnt.push_back(0);
                int nedges = edges.size();
                graphedges.emplace_back(lp.x, lp.y, rp.x, rp.y, i, -1);
                for (int j = 0; j < nedges - 1; j++) {
                    if ((edges[j].start == edges.back().start && edges[j].end == edges.back().end) ||
                        (edges[j].start == edges.back().end && edges[j].end == edges.back().start)) {
                        edgecnt[j]++;
                        edgecnt.back()++;
                        graphedges.back().site2 = j / 2;
                        graphedges[j / 2].site2 = i;
                    }
                }
            }
        }
    }

    bool check_circle_event(int i, int& left, int& right, double& x) {
        int n = sites.size();
        left = i - 1;
        right = i + 1;
        bool ileft = false, iright = false;
        double xleft, xright;
        while (left >= 0 && processed[left] == false) {
            double a = intersection_point(i, left);
            if (a > x) {
                xleft = a;
                ileft = true;
                break;
            }
            left--;
        }
        if (left < 0) {
            xleft = xmin;
        }
        while (right < n && processed[right] == false) {
            double a = intersection_point(i, right);
            if (a > x) {
                xright = a;
                iright = true;
                break;
            }
            right++;
        }
        if (right == n) {
            xright = xmax;
        }
        if (ileft && iright) {
            x = (xleft + xright) / 2;
            return true;
        }
        return false;
    }

    double intersection_point(int i, int j) {
        Site si = sites[i];
        Site sj = sites[j];
        double yi = si.p.y;
        double yj = sj.p.y;
        double xi = si.p.x;
        double xj = sj.p.x;
        double slopei = (xj - xi) / (yi - yj);
        double slopej = (xi - xj) / (yj - yi);
        double x = (si.p.x + sj.p.x + slopei*si.p.y - slopej*sj.p.y) / (slopei - slopej);
        return x;
    }

    void finish_edges() {
        int n = edges.size();
        for (int i = 0; i < n; i++) {
            if (edgecnt[i] == 1) {
                edges[i].end = vertices.size() > 0 ? vertices.back() : Point(NAN, NAN);
                vertices.push_back(Point(NAN, NAN));
            }
        }
    }

    void build_graph() {
        int nedges = edges.size();
        for (int i = 0; i < nedges; i++) {
            if (!isnan(edges[i].end.x)) {
                int s1 = edges[i].start.sitenum;
                int s2 = edges[i].end.sitenum;
                graphedges.emplace_back(edges[i].start.x, edges[i].start.y, edges[i].end.x, edges[i].end.y, s1, s2);
                neighbors[s1].push_back(graphedges.size() - 1);
                graphedges.emplace_back(edges[i].end.x, edges[i].end.y, edges[i].start.x, edges[i].start.y, s2, s1);
                neighbors[s2].push_back(graphedges.size() - 1);
            }
        }
    }

public:
    Voronoi() {}

    void compute() {
        initialize();
        insert_event(xmin - deltax, ymin - deltay);
        insert_event(xmin - deltax, ymax + deltay);
        insert_event(xmax + deltax, ymin - deltay);
        insert_event(xmax + deltax, ymax + deltay);
        process_event();
        finish_edges();
        build_graph();
    }

    vector<Point> get_vertices() {
        return vertices;
    }

    vector<GraphEdge> get_edges() {
        return graphedges;
    }

    vector<vector<int>> get_neighbors() {
        return neighbors;
    }
};

int main() {
    Voronoi v;
    v.compute();
    return 0;
}
``