[SPOJ CIRU]The Area of the Union of Circles

AC大神写过一篇圆面积并的文章,好像大家都是从他那边学来的(吐槽一下百度,自家空间的内容搜不到,但是谷歌就可以):求圆并的若干种算法,圆并扩展算法

首先,如果做过圆周长并([BZOJ1043][HAOI2008]下落的圆盘)的话,那做圆面积并就简单了。

根据AC大神的理论,圆的面积并的结果就是一大堆弓形+一个凸多边形(如果所有圆都交在一个块的话)。因此算面积并的话,就只要算出所有弓形的面积(1/2 * R * R * (α - sin(α))),还有所有有向线段的有向面积即可。所有的弓形其实就是每个圆没被覆盖的部分,这个比较好求。但是那个多边形呢?

AC大神提出了一种很厉害的方法:对于每个圆,枚举所有和它相交的圆,然后记下所有相交的区间,剩下的区间(若干条弧),计算弧所对应的弓形还有弧所对应的有向线段的有向面积。

ccir2

如图,当前圆为圆A,绿色和蓝色部分为当前圆没被覆盖的部分,而且一定是圆并中仅被覆盖一次的面积!实际上我们算的只有弧FG 弧GI向量FG 向量HI的面积,因为其他几条线段被包含在多边形EFGHI里面了,而且这个多边形的每条边都会被算到一次且仅一次。

然后有一种特殊情况:

ccir4

中间那块多边形明显是不算面积的,但是似乎是算了两次?如何解决BUG?

实际上这根本没有BUG,中间那块面积是负数!因为中间的有向线段方向是相反的!

所以这题就这样就可以做掉了(基本上就是圆周长并改一改)。代码也不长:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
#include <cmath>
#include <cstdio>
#include <vector>
#include <algorithm>
using std::pair;
using std::vector;

inline double sqr(const double x) { return x*x; }
inline double dist(const double x1, const double y1, const double x2, const double y2)
  { return std::sqrt(sqr(x1-x2) + sqr(y1-y2)); }
inline double angle(const double A, const double B, const double C)
  { return std::acos((sqr(A)+sqr(B)-sqr(C))/(2*A*B)); }
inline int sign(const double x)
{
  static const double EPS = 1e-8;
  if (x > EPS) return 1;
  return x < -EPS ? -1 : 0;
}
const int MAXN = 1000;
const double PI = std::acos(-1.0);
int n;
bool covered[MAXN];
double arc, pol, R[MAXN], A[MAXN], B[MAXN];
vector< pair<double, double> > seg, cover;
inline void getarea(const int i, const double lef, const double rig)
{
  arc += 0.5 * R[i] * R[i] * (rig-lef - std::sin(rig-lef));
  const double x1 = A[i] + R[i] * std::cos(lef), y1 = B[i] + R[i] * std::sin(lef);
  const double x2 = A[i] + R[i] * std::cos(rig), y2 = B[i] + R[i] * std::sin(rig);
  pol += x1*y2 - x2*y1;
}
int main()
{
  scanf("%d", &n);
  for (int i = 0; i < n; ++i)
  {
    scanf("%lf%lf%lf", A+i, B+i, R+i);
    for (int j = 0; j < i; ++j)
      if (!sign(R[i]-R[j]) && !sign(A[i]-A[j]) && !sign(B[i]-B[j]))
      {
        R[i] = .0;
        break;
      }
  }
  for (int i = 0; i < n; ++i)
    for (int j = 0; j < n; ++j)
      if (i != j && sign(R[j]-R[i]) >= 0 && sign(dist(A[i], B[i], A[j], B[j])-(R[j]-R[i])) <= 0)
      {
        covered[i] = true;
        break;
      }
  for (int i = 0; i < n; ++i) if (sign(R[i]) && !covered[i])
  {
    seg.clear();
    for (int j = 0; j < n; ++j) if (i != j && sign(R[i]) && !covered[i])
    {
      const double d = dist(A[i], B[i], A[j], B[j]);
      if (sign(d-(R[j]+R[i])) >= 0 || sign(d-std::abs(R[j]-R[i])) <= 0)
        continue;
      const double alpha = std::atan2(B[j]-B[i], A[j]-A[i]);
      const double beta = angle(R[i], d, R[j]);
      const pair<double, double> tmp(alpha-beta, alpha+beta);
      if (sign(tmp.first) <= 0 && sign(tmp.second) <= 0)
        seg.push_back(pair<double, double>(2*PI+tmp.first, 2*PI+tmp.second));
      else if (sign(tmp.first) < 0)
      {
        seg.push_back(pair<double, double>(2*PI+tmp.first, 2*PI));
        seg.push_back(pair<double, double>(0, tmp.second));
      }
      else
        seg.push_back(tmp);
    }

    std::sort(seg.begin(), seg.end());
    double rig = 0;
    for (vector< pair<double, double> >::iterator iter = seg.begin(); iter != seg.end(); ++iter)
    {
      if (sign(rig - iter->first) >= 0)
        rig = std::max(rig, iter->second);
      else
      {
        getarea(i, rig, iter->first);
        rig = iter->second;
      }
    }
    if (!sign(rig))
      arc += R[i]*R[i]*PI;
    else
      getarea(i, rig, 2*PI);
  }
  printf("%.3f", pol/2.0 + arc);
}

Comments