分割线算图形面积 圆和多边形面积并

冬令营听了hwd的课，感觉好像计算几何很简单的样子。

最近做了一道计算圆和多边形面积并的面积，用了扫面线的方法，感觉不太好（可能我没写习惯），写了\(7k\)的代码才A。QAQ

以后记得找月下柠檬树做做。

代码

#include <cstdio>
#include <iostream>
#include <vector>
#include <iomanip>
#include <cmath>
#include <assert.h>
using namespace std;

const int MAXN = 103;
const long double PI = acos(-1);
const long double EPS = 1e-15;

template <class T>
T sqr(const T &x) {
    return x * x;
}

int n, m;

struct Point {
    long double x, y;

    Point () {}

    Point (long double a, long double b):x(a), y(b) {}

    Point operator - (const Point &o) const {
        return Point(x - o.x, y - o.y);
    }

    Point operator + (const Point &o) const {
        return Point(x + o.x, y + o.y);
    }
    
    long double operator * (const Point &o) const {
        return x * o.y - y * o.x;
    }

    long double operator % (const Point &o) const {
        return x * o.x + y * o.y;
    }

    Point operator * (const long double &o) const {
        return Point(x * o, y * o);
    }

    Point operator / (const long double &o) const {
        return Point(x / o, y / o);
    }
};

struct Segment {
    Point first, second;

    Segment() {}

    Segment(Point x, Point y):first(x), second(y) {}
};

struct Circle {
    Point center;
    long double r;
};

Circle A[MAXN];
Point B[MAXN];

Point ret0, ret1;

long double dis(const Point &a, const Point &b) {
    return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y));
}

Point rotate(const Point &p, const long double &c, const long double &s) {
    return Point(p.x * c - p.y * s, p.x * s + p.y * c);
}

bool cir_cir(const Circle &a, const Circle &b) {
    long double len = dis(a.center, b.center);
    if (a.r + b.r < len) return false;
    if (len < fabs(a.r - b.r)) return false;
    long double COS = (sqr(a.r) + sqr(len) - sqr(b.r)) / (len * a.r * 2);
    long double SIN = sqrt(1. - sqr(COS));
    Point tmp = b.center - a.center;
    long double k = a.r / len;
    ret0 = a.center + rotate(tmp, COS, SIN) * k;
    ret1 = a.center + rotate(tmp, COS, -SIN) * k;
    return true;
}

Point line_line(const Segment &a, const Segment &b) {
    long double k = (a.second - a.first) * (b.first - a.first);
    k = k / (k - (a.second - a.first) * (b.second - a.first));
    return b.first + (b.second - b.first) * k;
}

bool cir_line(const Circle &a, const Segment &b) {
    Point v = b.second - b.first;
    swap(v.x, v.y);
    v.x = -v.x;
    v = a.center + v;
    Point A = line_line(Segment(v, a.center), b);
    long double d = dis(A - a.center, Point(0, 0));
    if (d > a.r) return false;
    long double len = sqrt(sqr(a.r) - sqr(d));
    long double norm = dis(b.second - b.first, Point(0, 0));
    ret0 = A + (b.second - b.first) * len / norm;
    ret1 = A - (b.second - b.first) * len / norm;
    return true;
}

long double cirSubArea(const Circle &a, const Point &p, const Point &q) {
    long double X = dis(p, q);
    long double angle = acos(1. - sqr(X) * .5 / sqr(a.r));
        //acos((sqr(a.r) * 2 - sqr(X)) * .5 / sqr(a.r));
    long double ret = angle * sqr(a.r) * .5;
    ret -= fabs((p - a.center) * (q - a.center) * .5);
    return ret;
}

struct Query {
    int belong;
    bool enter;
    Point left, right;
    long double key;

    bool operator < (const Query &o) const {
        return key < o.key;
    }
};

int main() {
    freopen("polygon.in", "r", stdin);
    freopen("polygon.out", "w", stdout);

    cin >> n >> m;
    vector<long double> key;

    for (int i = 0; i < n; i ++) {
        cin >> A[i].center.x >> A[i].center.y >> A[i].r;
        key.push_back(A[i].center.x - A[i].r);
        key.push_back(A[i].center.x + A[i].r);
        for (int j = 0; j < i; j ++) {
            if (cir_cir(A[i], A[j])) {
                key.push_back(ret0.x);
                key.push_back(ret1.x);
            }
        }
    }
    for (int i = 0; i < m; i ++) {
        cin >> B[i].x >> B[i].y;
        key.push_back(B[i].x);
    }
    B[m] = B[0];
    for (int i = 0; i < n; i ++) {
        for (int j = 0; j < m; j ++) {
            if (cir_line(A[i], Segment(B[j], B[j + 1]))) {
                if (ret0.x >= min(B[j].x, B[j + 1].x) - EPS &&
                        ret0.x <= max(B[j].x, B[j + 1].x) + EPS)
                    key.push_back(ret0.x);
                if (ret1.x >= min(B[j].x, B[j + 1].x) - EPS &&
                        ret1.x <= max(B[j].x, B[j + 1].x) + EPS)
                    key.push_back(ret1.x);
            }
        }
    }

    sort(key.begin(), key.end());
    vector<Query> que;

    long double ans = 0;

    for (int i = 0, _i = key.size(); i + 1 < _i; i ++) {
        if (i < _i && key[i] == key[i + 1]) continue;
        long double mid = (key[i] + key[i + 1]) * .5;
        Segment vecl = Segment(Point(key[i], 0), Point(key[i], 1));
        Segment vecr = Segment(Point(key[i + 1], 0), Point(key[i + 1], 1));
        Segment vecm = Segment(Point(mid, 0), Point(mid, 1));
        que.clear();

        for (int j = 0; j < n; j ++) {
            if (cir_line(A[j], vecm)) {
                if (ret0.y > ret1.y) swap(ret0, ret1);
                long double key0 = ret0.y, key1 = ret1.y;

                cir_line(A[j], vecl);
                Point a = ret0, b = ret1;
                if (a.y > b.y) swap(a, b);
                cir_line(A[j], vecr);
                if (ret0.y > ret1.y) swap(ret0, ret1);
                Query tmp;
                tmp.belong = j;

                tmp.left = a;
                tmp.right = ret0;
                tmp.enter = true;
                tmp.key = key0;
                que.push_back(tmp);
                tmp.left = b;
                tmp.right = ret1;
                tmp.enter = false;
                tmp.key = key1;
                que.push_back(tmp);
            }
        }
        bool in = false;
        for (int j = 0; j < m; j ++) {
            Point t = line_line(Segment(B[j], B[j + 1]), vecm);
            if (((B[j + 1] - B[j]) % (t - B[j]) > 0) == 
                    ((B[j + 1] - B[j]) % (t - B[j + 1]) > 0)) continue;

            in = true;
            Query tmp;
            tmp.belong = -1;
            tmp.left = line_line(vecl, Segment(B[j], B[j + 1]));
            tmp.right = line_line(vecr, Segment(B[j], B[j + 1]));
            tmp.enter = false;
            tmp.key = t.y;
            que.push_back(tmp);
        }
        if (in) {
            int t = que.size();
            if (que[t - 2].left.y >= que[t - 1].left.y &&
                    que[t - 2].right.y >= que[t - 1].right.y)
                que[t - 1].enter = true;
            else
                que[t - 2].enter = true;
        }

        sort(que.begin(), que.end());

        int cover = 0;
        Query lower;
        for (int j = 0, _j = que.size(); j < _j; j ++) {
            Query &cur = que[j];
            if (cur.enter) {
                if (cover == 0) {
                    lower = cur;
                    if (cur.belong != -1)
                        ans += cirSubArea(A[cur.belong], 
                                cur.left, cur.right);
                }
                cover ++;
            }
            else {
                cover --;
                if (cover == 0) {
                    ans += (cur.left.y + cur.right.y - lower.left.y - lower.right.y) *
                        (key[i + 1] - key[i]) * .5;
                    if (cur.belong != -1) 
                        ans += cirSubArea(A[cur.belong], cur.left, cur.right);
                }
            }
        }
    }

    //cout << ans << endl;
    printf("%.15lf\n", (double) ans);

    return 0;
}

转载于:https://www.cnblogs.com/wangck/p/4461534.html