SPOJ 8073 The area of the union of circles(计算几何の圆并)(CIRU)

Description

You are given N circles and expected to calculate the area of the union of the circles !

Input

The first line is one integer n indicates the number of the circles. (1 <= n <= 1000)

Then follows n lines every line has three integers

Xi Yi Ri

indicates the coordinate of the center of the circle, and the radius. (|Xi|. |Yi|  <= 1000, Ri <= 1000)

Note that in this problem Ri may be 0 and it just means one point !

Output

The total area that these N circles with 3 digits after decimal point

题目大意:求n个圆覆盖的总面积。

思路:参考http://hi.baidu.com/aekdycoin/item/b8ff6adc73c0e71dd78ed0d6

时间复杂度O(n^2*log(n))

代码(0.04S):

 #include <cstdio>
#include <algorithm>
#include <cstring>
#include <iostream>
#include <cmath>
#include <vector>
using namespace std;
typedef long long LL; const double PI = acos(-1.0);
const double EPS = 1e-; inline int sgn(double x) {
return (x > EPS) - (x < -EPS);
} struct Point {
double x, y;
Point() {}
Point(double x, double y): x(x), y(y) {}
void read() {
scanf("%lf%lf", &x, &y);
}
double angle() {
return atan2(y, x);
}
Point operator + (const Point &rhs) const {
return Point(x + rhs.x, y + rhs.y);
}
Point operator - (const Point &rhs) const {
return Point(x - rhs.x, y - rhs.y);
}
Point operator * (double t) const {
return Point(x * t, y * t);
}
Point operator / (double t) const {
return Point(x / t, y / t);
}
double length() const {
return sqrt(x * x + y * y);
}
Point unit() const {
double l = length();
return Point(x / l, y / l);
}
}; double cross(const Point &a, const Point &b) {
return a.x * b.y - a.y * b.x;
} double dist(const Point &p1, const Point &p2) {
return (p1 - p2).length();
} Point rotate(const Point &p, double angle, const Point &o = Point(, )) {
Point t = p - o;
double x = t.x * cos(angle) - t.y * sin(angle);
double y = t.y * cos(angle) + t.x * sin(angle);
return Point(x, y) + o;
} struct Region {
double st, ed;
Region() {}
Region(double st, double ed): st(st), ed(ed) {}
bool operator < (const Region &rhs) const {
if(sgn(st - rhs.st)) return st < rhs.st;
return ed < rhs.ed;
}
}; struct Circle {
Point c;
double r;
vector<Region> reg;
Circle() {}
Circle(Point c, double r): c(c), r(r) {}
void read() {
c.read();
scanf("%lf", &r);
}
void add(const Region &r) {
reg.push_back(r);
}
bool contain(const Circle &cir) const {
return sgn(dist(cir.c, c) + cir.r - r) <= ;
}
bool intersect(const Circle &cir) const {
return sgn(dist(cir.c, c) - cir.r - r) < ;
}
}; double sqr(double x) {
return x * x;
} void intersection(const Circle &cir1, const Circle &cir2, Point &p1, Point &p2) {
double l = dist(cir1.c, cir2.c);
double d = (sqr(l) - sqr(cir2.r) + sqr(cir1.r)) / ( * l);
double d2 = sqrt(sqr(cir1.r) - sqr(d));
Point mid = cir1.c + (cir2.c - cir1.c).unit() * d;
Point v = rotate(cir2.c - cir1.c, PI / ).unit() * d2;
p1 = mid + v, p2 = mid - v;
} Point calc(const Circle &cir, double angle) {
Point p = Point(cir.c.x + cir.r, cir.c.y);
return rotate(p, angle, cir.c);
} const int MAXN = ; Circle cir[MAXN];
bool del[MAXN];
int n; double solve() {
double ans = ;
for(int i = ; i < n; ++i) {
for(int j = ; j < n; ++j) if(!del[j]) {
if(i == j) continue;
if(cir[j].contain(cir[i])) {
del[i] = true;
break;
}
}
}
for(int i = ; i < n; ++i) if(!del[i]) {
Circle &mc = cir[i];
Point p1, p2;
bool flag = false;
for(int j = ; j < n; ++j) if(!del[j]) {
if(i == j) continue;
if(!mc.intersect(cir[j])) continue;
flag = true;
intersection(mc, cir[j], p1, p2);
double rs = (p2 - mc.c).angle(), rt = (p1 - mc.c).angle();
if(sgn(rs) < ) rs += * PI;
if(sgn(rt) < ) rt += * PI;
if(sgn(rs - rt) > ) mc.add(Region(rs, PI * )), mc.add(Region(, rt));
else mc.add(Region(rs, rt));
}
if(!flag) {
ans += PI * sqr(mc.r);
continue;
}
sort(mc.reg.begin(), mc.reg.end());
int cnt = ;
for(int j = ; j < int(mc.reg.size()); ++j) {
if(sgn(mc.reg[cnt - ].ed - mc.reg[j].st) >= ) {
mc.reg[cnt - ].ed = max(mc.reg[cnt - ].ed, mc.reg[j].ed);
} else mc.reg[cnt++] = mc.reg[j];
}
mc.add(Region());
mc.reg[cnt] = mc.reg[];
for(int j = ; j < cnt; ++j) {
p1 = calc(mc, mc.reg[j].ed);
p2 = calc(mc, mc.reg[j + ].st);
ans += cross(p1, p2) / ;
double angle = mc.reg[j + ].st - mc.reg[j].ed;
if(sgn(angle) < ) angle += * PI;
ans += 0.5 * sqr(mc.r) * (angle - sin(angle));
}
}
return ans;
} int main() {
scanf("%d", &n);
for(int i = ; i < n; ++i) cir[i].read();
printf("%.3f\n", solve() + EPS);
}

以下转自:http://hi.baidu.com/aekdycoin/item/b8ff6adc73c0e71dd78ed0d6

【求圆并的若干种算法,圆并扩展算法】

【问题求解】
给定N 个圆形,求出其并集.

【算法分析】

PS.以下算法基于正方向为逆时针

SPOJ 8073 The area of the union of circles(计算几何の圆并)(CIRU)

考虑上图中的蓝色圆,绿色的圆和蓝色的圆交于 A,B 2个交点 ,我们在逆时针系下考虑,那么 可以知道 对于蓝色的圆,它对应的某个 角度区间被覆盖了

假设 区间为 [A, B], 并且角度是按照 圆心到交点的 向量的 极角来定义 (为了方便,我一般都把角度转化到 [0,2pi]区间内) 那么可以知道在这种 标识情况下,可能存在以下情况

SPOJ 8073 The area of the union of circles(计算几何の圆并)(CIRU)

这种区间的跨度如何解决呢?实际上十分简单,只需要把[A,B] 拆成 [A, 2PI], [0,B]即可,也就是所谓的添加虚拟点

SPOJ 8073 The area of the union of circles(计算几何の圆并)(CIRU)

下面介绍一下 对于我们当前所求任务的实际运用( 利用上述做法)

首先 对于所给的N个圆,我们可以进行去冗杂,表现为:
(1) 去掉被包含(内含 or 内切)的小圆 ()
(2) 去掉相同的圆

枚举一个圆,并对于剩下的圆和它求交点,对于所求的的交点,可以得到一个角度区间 [A,B], 当然区间如果跨越了(例如 [1.5PI, 0.5PI],注意这里是有方向的) 2PI那么需要拆 区间

可以知道,最后区间的并集必然是最后 所有圆和当前圆的交集的一个边界!

于是我们得到互补区间(所谓互补区间就是[0,2PI]内除去区间并的区间,可能有多个)

假设我们先枚举了橙色的圆,那么得到了许多角度区间,可以知道绿色的和蓝色的角度区间是“未被覆盖的”,对于未被覆盖的

圆弧染色!

而对于其他圆,我们也做同样的步骤, 同时把他们未被覆盖的角度区间的圆弧标记为黑色阴影

于是最终的结果就如下图 (染色只考虑圆弧)

SPOJ 8073 The area of the union of circles(计算几何の圆并)(CIRU)

通过观察不难发现,圆的并是许多的圆弧+ 许多多边形的面积之和(注意这里为简单多边形,并且面积有正负之别!)

于是我们累加互补区间的圆弧面积到答案,并把互补区间确定的弦的有向面积累加到答案

(例如在上面的例子中,我们在扫描到橙色圆的这一步只需要累加蓝色和绿色的圆弧面积 以及 蓝色和绿色的有向面积,注意这里蓝色和绿色的边必然是最后那个多边形的边!)

这里涉及到一个问题,就是:
圆弧可能大于一半的圆,例如上图中最大的圆,当然如果你推出了公式,那么实际上很容易发现无论圆弧长啥样都能算出正确的答案!
半径为R的圆中,弧度区间为K的圆弧的面积为 : 0.5 * r * r * (K - sin(K))

之后还是有一些疑惑,如果存在"洞"?

美妙之处在于此,由于面积有方向(正负),所以洞会被自然抵消!(可以类比计算简单多边形面积)

枚举 圆,并和其他圆求交,求得区间,排序
这里是O(n^2 logn)
至于区间扫描和累加的复杂度则是O(n)
于是复杂度
O(n * (nlogn + n) )

O(n^2 logn)

代码已实现,大约140行,和圆的交的思路十分类似,圆的交的思路请参考:

http://hi.baidu.com/aekdycoin/blog/item/12267a4e9476153bafc3abbd.html

圆并的题目:

http://www.spoj.pl/problems/CIRU/

https://www.spoj.pl/problems/VCIRCLES/

【圆并的数值积分】

// 先贴代码,稍后补上

 #include<iostream>
#include<stdio.h>
#include<string.h>
#include<algorithm>
#include<cmath>
using namespace std;
const int maxn = ;
typedef double db;
const db EPS = 1e-;
typedef pair<db, db> PDD;
int x[ maxn ], y[ maxn ], r[ maxn ];
int nx[ maxn ], ny[ maxn ], nr[ maxn ];
int xl[ maxn ], xr[ maxn ]; int s[ maxn ];
inline bool cmp( int a, int b) {
if( x[ a ] - r [ a ] == x[ b ] - r [ b ] ) return x[ a ] + r[ a ] < x[ b ] + r [ b ];
return x[ a ] - r[ a ] < x[ b ] - r [ b ];
}
inline bool cmp0(int a, int b){return r[ a ] > r [ b ];}
int n;
int L, R;
PDD se[ maxn ];
inline db f( db v){
int sz = , i, j ;
db ret = 0.0;
for(i = L; i < R; ++ i){
if( v <= xl[ i ] || v >= xr[ i ] ) continue;
j = s[ i ];
db d = sqrt(r[ j ]- (v - x [ j ]) * (v - x[ j ]));
se[ sz ].first = y[ j ] - d;
se[ sz ].second = y[ j ] + d;
++ sz;
}
sort( se, se + sz);
for(i = ; i < sz; ++ i){
db nowl , nowr;
nowl = se[ i ].first;
nowr = se[ i ].second;
for( j = i + ; j < sz; ++ j) if(se[ j ].first > nowr) break;
else nowr = max( nowr, se[ j ].second);
ret += nowr - nowl;
i = j - ;
}
return ret;
}
#define fs(x) ((x) < 0 ? (-(x)) : (x))
inline db rsimp( db l,db m, db r, db sl, db sm, db sr,db tot){
db m1 = (l + m) * 0.5, m2 = (m + r) * 0.5;
db s0 = f( m1), s2 = f( m2);
db gl = (sl + sm + s0 + s0 + s0 + s0)*(m-l), gr = (sm + sr + s2 + s2 + s2 + s2)*(r-m);
if( fs(gl + gr - tot) < EPS) return gl + gr;
return rsimp( l, m1, m, sl, s0, sm, gl) + rsimp( m, m2,r, sm, s2, sr, gr);
} bool get(){
if( != scanf("%d", &n)) return ;
int i, j = , k;
for(i = ; i < n; ++ i) scanf("%d%d%d", x + i, y + i, r + i), s[ i ] = i;
sort( s, s + n, cmp0);
for(i = ; i < n; ++ i){
for(k = ; k < j; ++ k)
if( (nx [ k ] - x [s[i]]) * (nx [ k ] - x [s[i]]) + (ny [ k ] - y [s[i]]) *(ny [ k ] - y [s[i]])
<= (nr[ k ] - r[ s[ i ] ]) * (nr[ k ] - r[ s[ i ] ]) ) break;
if( k == j) {
nx[ j ] = x[ s[ i ] ];
ny[ j ] = y[ s[ i ] ];
nr[ j ] = r[ s[ i ] ];
s[ j ] = j;
j ++;
}
}
n = j;
for(i = ; i < n; ++ i) x[ i ] = nx[ i ], y[ i ] = ny[ i ], r[ i ] = nr[ i ];
return ;
} void work(){
sort( s, s + n, cmp) ;
db lo, hi, ans = 0.0;
int i, j;
for(i = ; i < n; ++ i) xl[ i ] = x[ s[ i ] ] - r[ s[ i ] ], xr[ i ] = x[ s[ i ] ] + r[ s[ i ] ], r[ s[i] ] *= r[ s[i] ];
for(i = ; i < n; ++ i){
int ilo, ihi;
ilo = xl[ i ];
ihi = xr[ i ];
for(j = i + ; j < n; ++ j) {
if( xl[ j ] > ihi ) break;
ihi = max( ihi, xr[ j ]);
}
db lo = ilo;
db hi = ihi;
L = i;
R = j;
db mid = (lo + hi) * 0.5;
db sl = f(lo), sm = f(mid), sr = f(hi);
db tot = sl + sr + sm + sm + sm + sm;
ans += rsimp( lo, mid , hi, sl, sm , sr, tot );
i = j - ;
}
printf("%.3f\n", ans / 6.0);
} int main(){
while( get() ) work();
return ;
}

【圆并的扩展算法】

https://www.spoj.pl/problems/CIRUT/

【题目大意】

给出N个不同的圆,求出被覆盖恰好K次的面积,并顺序输出

SPOJ 8073 The area of the union of circles(计算几何の圆并)(CIRU)

【思路提示】

参考圆并和圆交的思路, 就是把问题转化为区间来考虑

同时请看下面的图:

SPOJ 8073 The area of the union of circles(计算几何の圆并)(CIRU)

上面的图,是对于被覆盖恰好2次的区域的边界连接起来得到的图形,我们可以发现什么呢?

(就不剧透了~~~~)


上一篇:IE8中jQuery.load()加载页面不显示的原因


下一篇:中频IF