【SPOJ CUIR】圆的面积并

【参考资料】

https://www.cnblogs.com/heisenberg-/p/6740654.html

 

【前置知识】

 

 

 

---

【题目见标题】

【题目大意】给 n 个圆，给其圆心何半径，求这 n 个圆覆盖的面积

【算法分析】

两种方法，一种是几何法，一种是自适应辛普森

 

之前在看自适应辛普森，结果卡时间卡的很不舒服，看底下几何法过的都闪电快如飞，于是决定用几何法

 

【几何法】

 通过计算圆和圆之间的交点，通过交点将圆的面积分为没有覆盖的圆和多次覆盖的多边形

然后将交点加入队列，极角排序后计算

通过对面积的计算，将结果变为仅仅计算一次

 

 

 

 

【代码】

 

#include<cstdio>
#include<cmath>
#include<vector>
#include<algorithm>
#include<map>
#include<iostream>
#include<complex>
#include<queue>
#include<stack>
using namespace std;
const double eps  = 1e-9;
const double PI   = acos(-1.0);
const int    MAXN = 1010000;
int sgn(double x)
{
    if (x < -eps)    return -1;
    if (x > eps)    return 1;
    return 0;
}
struct V
{
    double x, y;
    double ang;
    V(double X = 0, double Y = 0)
    {
        //初始化
        x = X, y = Y;
    }
    bool operator ==(const V &b)
    {
        return sgn(x - b.x) && sgn(y - b.y);
    }

};
V operator +(V a, V b) { return V(a.x + b.x, a.y + b.y); }
V operator -(V a, V b) { return V(a.x - b.x, a.y - b.y); }
V operator *(V a, double b) { return V(a.x *b, a.y*b); }
V operator /(V a, double b) { return V(a.x / b, a.y / b); }
//叉积
double cross(V a, V b)
{
    return a.x*b.y - a.y*b.x;
}
typedef V P;
//点积
struct node
{
    P p;
    double ang;
    int id;
    node(P a, double b, int c) :p(a), ang(b), id(c) {}
    bool operator<(node a)
    {
        return ang < a.ang;
    }
};
double dot(V a, V b)
{
    return a.x*b.x + a.y*b.y;
}
struct L
{
    P s, t;
    L() {}
    L(P a, P b):s(a), t(b) {}
};
typedef L  S;
struct C
{
    P p;
    double r;
    bool input()
    {
        return scanf("%lf%lf", &p.x, &p.y) != EOF;
    }
    bool operator<(const C &o)
    {
        if (sgn(r - o.r)     != 0)    return sgn(r - o.r) == -1;
        if (sgn(p.x - o.p.x) != 0)    return sgn(p.x - o.r) == -1;
        return sgn(p.y - o.p.y) == -1;
    }
    bool operator==(const C &o)
    {
        return sgn(r - o.r) == 0 && sgn(p.x - o.p.x)==0 && sgn(p.y - o.p.y)==0;
    }
};
double len(P a)
{
    return sqrt(a.x*a.x + a.y*a.y);
}
C c[MAXN],arr[MAXN];
int N;
P rotate(P a, double cost, double sint)
{
    double x = a.x, y = a.y;
    return P(x*cost - y * sint, x*sint + y * cost);
}
pair<P, P> crossP(P a, double r1, P b, double r2)
{
    double d = len(a - b);
    //计算两圆心的距离
    double cost = (r1*r1 + d * d - r2 * r2) / (2 * r1 *d);
    //利用相似求出
    double sint=(sqrt(1. - cost * cost));
    P v = (b - a) / len(b - a) * r1;
    return make_pair(a + rotate(v, cost,-sint), a + rotate(v, cost, sint));
    
}
pair<P, P> crossP(C a, C b)
{
    return crossP(a.p, a.r, b.p, b.r);
    //方便填写函数的重载
}
double abs_t(double x)
{
    if (x < -eps) return -x;
    return x;
}
double arg(P a)
{
    return arg(complex<double>(a.x, a.y));
}
vector<node> vec;
double solve()
{
    int n = 0;
    sort(c+1, c + 1 + N);
    N = unique(c + 1, c + 1 + N) - (c + 1);
    for (int i = N; i >= 1; i--)
    {
        bool flag = true;
        for (int j = i + 1; j <= N; j++)
        {
            double d = len((c[i].p - c[j].p));
            if (sgn(d - abs_t(c[i].r - c[j].r)) <= 0)
            {
                flag = false; break;
            }
        }
        if (flag) arr[++n] = c[i];
        //去除被包含的圆
    }
    double ans = 0;
    for (int i = 1; i <= n; i++)
    {
        vec.clear();
        P bound = arr[i].p + P(-arr[i].r, 0);
        vec.push_back(node(bound, -PI, 0));
        vec.push_back(node(bound, PI , 0));
        for (int j = 1; j <= n; j++)
        {
            if (i == j)    continue;
            double d = len(arr[i].p - arr[j].p);
            if (sgn(d - (arr[i].r + arr[j].r)) < 0)
            {
                pair<P, P> ret = crossP(arr[i], arr[j]);
                //求取到两个圆的交点
                double x = arg(ret.first  - arr[i].p);
                //arg函数，求这个交点相对圆心的角度
                double y = arg(ret.second - arr[i].p);
                if (sgn(x - y) > 0)
                {
                    vec.push_back(node(ret.first , x , 1 ));
                    vec.push_back(node(bound     , PI, -1));
                    vec.push_back(node(bound     , -PI, 1 ));
                    vec.push_back(node(ret.second, y , -1));
                }
                else
                {
                    vec.push_back(node(ret.first , x, 1 ));
                    vec.push_back(node(ret.second, y, -1));
                }
            }
        }
        sort(vec.begin(), vec.end());
        int sum = vec[0].id;
        for (int j = 1; j <(int)vec.size(); j++)
        {
            if (sum==0) 
            {
                ans += cross(vec[j - 1].p, vec[j].p) / 2.;
                double x = vec[j - 1].ang;
                double y = vec[j].ang;
                double area = arr[i].r*arr[i].r*(y - x) / 2.;
                P v1 = vec[j - 1].p - arr[i].p;
                P v2 = vec[j].p - arr[i].p;
                area -= cross(v1, v2) / 2.;
                ans += area;
            }
            sum += vec[j].id;
        }
    }
    return ans;
}
int main()
{
    int m;
    scanf("%d", &m);
    N = 0;
    for (int i = 1; i <= m; i++)
    {
        C tmp;
        tmp.input();
        scanf("%lf",&tmp.r);
        if (sgn(tmp.r) != 0)
            c[++N] = tmp;
    }
    printf("%.3f\n",solve());
    return 0;
}

 

 

 

【自适应辛普森】

 利用积分的性质，计算每个圆弧