整理了一下计算几何的模板
后续会继续更新。。。。。
const double eps = 1e-8;
const double PI = acos(-1.0);
int sgn(double x){
if(fabs(x) < eps)return 0;
if(x < 0)return -1;
else return 1;
}
//结构体定义
struct Point {
double x,y;
Point(){}
Point(double xx,double yy){x=xx;y=yy;}
Point operator-(Point b){return Point(x-b.x,y-b.y);}
Point operator+(Point b){return Point(x+b.x,y+b.y);}
double operator^(Point b){return x*b.y-y*b.x;}
double operator*(Point b){return x*b.x-y*b.y;}
//绕原点逆时针转B度
void transXY(double B){
double tx = x,ty = y;
x = tx*cos(B) - ty*sin(B);
y = tx*sin(B) + ty*cos(B);
}
}p[maxn];
struct Line {
Point s,e;
Line(){}
Line(Point ss,Point ee){s=ss;e=ee;}
//两直线求交点
pair<int,Point> operator &(const Line &b)const{
Point res = s;
if(sgn((s-e)^(b.s-b.e)) == 0){
if(sgn((s-b.e)^(b.s-b.e)) == 0)
return make_pair(0,res);//重合
else return make_pair(1,res);//平行
}
double t = ((s-b.s)^(b.s-b.e))/((s-e)^(b.s-b.e));
res.x += (e.x-s.x)*t;
res.y += (e.y-s.y)*t;
return make_pair(2,res);
}
};
double dist(Point a,Point b){
return sqrt((a-b)*(a-b));
}
bool inter(Line l1,Line l2){
return
max(l1.s.x,l1.e.x) >= min(l2.s.x,l2.e.x) &&
max(l2.s.x,l2.e.x) >= min(l1.s.x,l1.e.x) &&
max(l1.s.y,l1.e.y) >= min(l2.s.y,l2.e.y) &&
max(l2.s.y,l2.e.y) >= min(l1.s.y,l1.e.y) &&
sgn((l2.s-l1.e)^(l1.s-l1.e))*sgn((l2.e-l1.e)^(l1.s-l1.e)) <= 0 &&
sgn((l1.s-l2.e)^(l2.s-l2.e))*sgn((l1.e-l2.e)^(l2.s-l2.e)) <= 0;
}
//判断直线和线段相交
bool Seg_inter_line(Line l1,Line l2) {//判断直线l1和线段l2是否相交
return sgn((l2.s-l1.e)^(l1.s-l1.e))*sgn((l2.e-l1.e)^(l1.s-l1.e)) <= 0;
}
//点到直线距离
//返回为result,是点到直线最近的点
Point PointToLine(Point P,Line L){
Point result;
double t = ((P-L.s)*(L.e-L.s))/((L.e-L.s)*(L.e-L.s));
result.x = L.s.x + (L.e.x-L.s.x)*t;
result.y = L.s.y + (L.e.y-L.s.y)*t;
return result;
}
//点到线段的距离
//返回点到线段最近的点
Point NearestPointToLineSeg(Point P,Line L){
Point result;
double t = ((P-L.s)*(L.e-L.s))/((L.e-L.s)*(L.e-L.s));
if(t >= 0 && t <= 1){
result.x = L.s.x + (L.e.x - L.s.x)*t;
result.y = L.s.y + (L.e.y - L.s.y)*t;
}
else{
if(dist(P,L.s) < dist(P,L.e))
result = L.s;
else result = L.e;
}
return result;
}
//计算多边形面积
//点的编号从0~n-1
double CalcArea(Point p[],int n){
double res = 0;
for(int i = 0;i < n;i++)
res += (p[i]^p[(i+1)%n])/2;
return fabs(res);
}
//*判断点在线段上
bool OnSeg(Point P,Line L){
return
sgn((L.s-P)^(L.e-P)) == 0 &&
sgn((P.x - L.s.x) * (P.x - L.e.x)) <= 0 &&
sgn((P.y - L.s.y) * (P.y - L.e.y)) <= 0;
}
//*判断点在凸多边形内
//点形成一个凸包，而且按逆时针排序（如果是顺时针把里面的<0改为>0）
//点的编号:0~n-1
//返回值：
//-1:点在凸多边形外
//0:点在凸多边形边界上
//1:点在凸多边形内
int inConvexPoly(Point a,Point p[],int n){
for(int i = 0;i < n;i++){
if(sgn((p[i]-a)^(p[(i+1)%n]-a)) < 0)return -1;
else if(OnSeg(a,Line(p[i],p[(i+1)%n])))return 0;
}
return 1;
}
//*判断点在任意多边形内
//射线法，poly[]的顶点数要大于等于3,点的编号0~n-1
//返回值
//-1:点在凸多边形外
//0:点在凸多边形边界上
//1:点在凸多边形内
int inPoly(Point p,Point poly[],int n){
int cnt;
Line ray,side;
cnt = 0;
ray.s = p;
ray.e.y = p.y;
ray.e.x = -100000000000.0;//-INF,注意取值防止越界
for(int i = 0;i < n;i++){
side.s = poly[i];
side.e = poly[(i+1)%n];
if(OnSeg(p,side))return 0;
//如果平行轴则不考虑
if(sgn(side.s.y - side.e.y) == 0)
continue;
if(OnSeg(side.s,ray))
{
if(sgn(side.s.y - side.e.y) > 0)cnt++;
}
else if(OnSeg(side.e,ray))
{
if(sgn(side.e.y - side.s.y) > 0)cnt++;
}
else if(inter(ray,side))
cnt++;
}
if(cnt % 2 == 1) return 1;
else return -1;
}
//判断凸多边形
//允许共线边
//点可以是顺时针给出也可以是逆时针给出
//点的编号1~n-1
bool isconvex(Point poly[],int n){
bool s[3];
memset(s,false,sizeof(s));
for(int i = 0;i < n;i++)
{
s[sgn( (poly[(i+1)%n]-poly[i])^(poly[(i+2)%n]-poly[i]) )+1] = true;
if(s[0] && s[2])return false;
}
return true;
}
/*
* 求凸包，Graham算法
* 点的编号0~n-1
* 返回凸包结果ch[0~top-1]为凸包的编号
* poj2007 凸包点排序逆时针输出
*/
#include <stdio.h>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
#define sqr(a) ((a) * (a))
#define dis(a, b) sqrt(sqr(a.x - b.x) + sqr(a.y - b.y))
const int MAXN = 110;
const double PI = acos(-1.0);
struct Point {
int x;
int y;
Point(double a = 0, double b = 0) : x(a), y(b) {}
friend bool operator < (const Point &l, const Point &r) {
return l.y < r.y || (l.y == r.y && l.x < r.x);
}
} p[MAXN], ch[MAXN];
double mult(Point a, Point b, Point o) {
return (a.x - o.x) * (b.y - o.y) >= (b.x - o.x) * (a.y - o.y);
}
int Graham(Point p[], int n, Point res[]) {
int top = 1;
sort(p, p + n);
if (n == 0) return 0;
res[0] = p[0];
if (n == 1) return 0;
res[1] = p[1];
if (n == 2) return 0;
res[2] = p[2];
int i;
for (i = 2; i < n; i++) {
while (top && (mult(p[i], res[top], res[top - 1])))
top--;
res[++top] = p[i];
}
int len = top;
res[++top] = p[n - 2];
for (i = n - 3; i >= 0; i--) {
while (top != len && (mult(p[i], res[top], res[top - 1])))
top--;
res[++top] = p[i];
}
return top;
}
int n;
int main() {
while (scanf("%d%d", &p[n].x, &p[n].y) != EOF)
n++;
n = Graham(p, n, ch);
int t,i;
for (i = 0; i < n; i++)
if (ch[i].x == 0 && ch[i].y == 0) {
t = i;
break;
}
for (i = t; i < n; i++)
printf("(%d,%d)\n", ch[i].x, ch[i].y);
for (i = 0; i < t; i++)
printf("(%d,%d)\n", ch[i].x, ch[i].y);
return 0;
}
/*
平面最近点对 ZOJ 2107
*/
#include<stdio.h>
#include<math.h>
#include<algorithm>
using namespace std;
const double eps=1e-8;
struct Point {
double x,y;
}p[100005],tmp[100005];
int n;
double dist(Point a,Point b){
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
int cmpxy(Point a,Point b){
if(a.x!=b.x) return a.x<b.x;
return a.y<b.y;
}
int cmpy(Point a,Point b){
return a.y<b.y;
}
double Closest_Pair(int left,int right){
double d=1e20;
if(left==right) return d;
if(left+1==right) return dist(p[left],p[right]);
int mid=(left+right)>>1;
double d1=Closest_Pair(left,mid);
double d2=Closest_Pair(mid+1,right);
d=min(d1,d2);
int k=0;
for(int i=left;i<=right;i++){
if(fabs(p[mid].x-p[i].x)<d+eps) tmp[k++]=p[i];
}
sort(tmp,tmp+k,cmpy);
for(int i=0;i<k;i++){
for(int j=i+1;j<k&&tmp[j].y-tmp[i].y<d+eps;j++){
d=min(d,dist(tmp[i],tmp[j]));
}
}
return d;
}
int main(){
int n;
while(scanf("%d",&n),n){
for(int i=0;i<n;i++) scanf("%lf%lf",&p[i].x,&p[i].y);
sort(p,p+n,cmpxy);
printf("%.2lf\n",Closest_Pair(0,n-1)/2);
}
}
/*
旋转卡壳
求平面最远点对 Beauty Contest POJ 2187
*/
#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <cmath>
#include<algorithm>
#define N 50005
using namespace std;
const double eps=1e-5;
struct Point {
int x,y;
Point (){}
Point (int xx,int yy) {x=xx;y=yy;}
Point operator -(const Point b)const {return Point(x-b.x,y-b.y);}
int operator ^(const Point b)const {return x*b.y-y*b.x;}
};
int n,top;
Point p[N],stack[N];
int dist(Point a,Point b){
return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);
}
int cmp(Point a,Point b) {//逆时针
int tmp=(a-p[0])^(b-p[0]);
if(tmp == 0 ) return dist(a,p[0]) - dist(b,p[0]) <= 0;
else return tmp>0;
}
void Graham(){
int k=0;
for(int i=0;i<n;i++){
if(p[i].y<p[k].y || (p[i].y==p[k].y && p[i].x<p[k].x)) k=i;
}
swap(p[0],p[k]);
sort(p+1,p+n,cmp);
if(n<=2){
for(int i=0;i<n;i++) stack[i]=p[i];
top=n;
return;
}
stack[0]=p[0]; stack[1]=p[1];
top=2;
for(int i=2;i<n;i++){
while(top>=2 && ((stack[top-1]-stack[top-2])^(p[i]-stack[top-2]))<=0) top--;
stack[top++]=p[i];
}
}
int main(){
while(scanf("%d",&n)!=EOF){
for(int i=0;i<n;i++) scanf("%d%d",&p[i].x,&p[i].y);
Graham();
for(int i=0;i<top;i++) {
p[i]=stack[i];
}
p[top]=p[0];
int ans=0;
Point a,b;
for(int i=0,j=1;i<top;i++){
while(((p[i+1]-p[i])^(p[j]-p[i]))<((p[i+1]-p[i])^(p[j+1]-p[i]))) j=(j+1)%top;
if(dist(p[i],p[j])>=ans) ans=dist(p[i],p[j]);
if(dist(p[i+1],p[j+1])>=ans) ans=dist(p[i+1],p[j+1]);//处理两条边平行的特殊情况
}
printf("%d\n",ans);
}
return 0;
}
//旋转卡壳计算平面点集最大三角形面积
int rotating_calipers(Point p[],int n){
int ans = 0;
Point v;
for(int i = 0;i < n;i++){
int j = (i+1)%n;
int k = (j+1)%n;
while(j != i && k != i){
ans = max(ans,abs((p[j]-p[i])^(p[k]-p[i])));
while( ((p[i]-p[j])^(p[(k+1)%n]-p[k])) < 0 )
k = (k+1)%n;
j = (j+1)%n;
}
}
return ans;
}
Point p[MAXN];
int main(){
int n;
while(scanf("%d",&n) == 1){
if(n == -1)break;
for(int i = 0;i < n;i++)list[i].input();
Graham(n);
for(int i = 0;i < top;i++)p[i] = list[Stack[i]];
printf("%.2f\n",(double)rotating_calipers(p,top)/2);
}
return 0;
}
/*
最近点对和最远点对产生在有向平行切线上，将一条边平行过去，然后判断
求两凸包最近点对，最远点对类似
*/
#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <cmath>
#include<algorithm>
using namespace std;
const double eps=1e-8;
struct Point {
double x,y;
Point (){}
Point (double xx,double yy) {x=xx;y=yy;}
Point operator -(const Point b)const {return Point(x-b.x,y-b.y);}
double operator ^(const Point b)const {return x*b.y-y*b.x;}
double operator *(const Point &b)const{
return x*b.x + y*b.y;
}
};
struct Line {
Point s,e;
double angle;
Line(){}
Line(Point ss,Point ee) {s=ss;e=ee;}
};
Point p[10005],q[10005],list[10005];
int n,m,top,miny,maxy;
double dist(Point a,Point b) {
return sqrt((b.x-a.x)*(b.x-a.x)+(b.y-a.y)*(b.y-a.y));
}
//求点到线段距离，不是点到直线距离
double pointtoseg(Point P,Point b,Point c){
Point result;
Line L=Line(b,c);
double t = ((P-L.s)*(L.e-L.s))/((L.e-L.s)*(L.e-L.s));
if(t >=0 && t <= 1)
{
result.x = L.s.x + (L.e.x - L.s.x)*t;
result.y = L.s.y + (L.e.y - L.s.y)*t;
}
else
{
if(dist(P,L.s) < dist(P,L.e))
result = L.s;
else result = L.e;
}
return dist(result,P);
}
double segtoseg(Point a,Point b,Point c,Point d){
double ans1=min(pointtoseg(a,c,d),pointtoseg(b,c,d));
double ans2=min(pointtoseg(c,a,b),pointtoseg(d,a,b));
return min(ans1,ans2);
}
int cmp(Point a,Point b) {//逆时针
double tmp=(a-list[0])^(b-list[0]);
if(tmp == 0) return dist(a,list[0]) - dist(b,list[0]) <-eps;
else return tmp>0;
}
void Graham(Point* p,Point* stack,int n){
int k=0;
for(int i=0;i<n;i++){
if(p[i].y<p[k].y || (p[i].y==p[k].y && p[i].x<p[k].x)) k=i;
}
swap(p[0],p[k]);
sort(p+1,p+n,cmp);
if(n<=2){
for(int i=0;i<n;i++) stack[i]=p[i];
top=n;
return;
}
stack[0]=p[0]; stack[1]=p[1];
top=2;
for(int i=2;i<n;i++){
while(top>=2 && ((stack[top-1]-stack[top-2])^(p[i]-stack[top-2]))<eps) top--;
stack[top++]=p[i];
}
}
double Get_angle(Point a1,Point a2,Point b1,Point b2)
{
return (a2-a1)^(b1-b2);
}
double rotating(Point* p,int n,Point* q,int m){
miny=0,maxy=0;
for(int i=0;i<n;i++) {
if(p[i].y<p[miny].y) miny=i;
}
for(int i=0;i<m;i++) {
if(q[i].y>q[maxy].y) maxy=i;
}
double ans=dist(p[miny],q[maxy]);
for(int i=0;i<n;i++){
double tmp;
while((tmp=Get_angle(p[miny],p[miny+1],q[maxy],q[maxy+1]))<-eps) maxy=(maxy+1)%m;
if(fabs(tmp)<eps) ans=min(ans,segtoseg(p[miny],p[(miny+1)%n],q[maxy],q[(maxy+1)%m]));
else ans=min(ans,pointtoseg(q[maxy],p[miny],p[(miny+1)%n]));
miny=(miny+1)%n;
}
return ans;
}
int main() {
while(scanf("%d%d",&n,&m)) {
if(n==0 && m==0) break;
//顺时针排列的点
miny=0,maxy=0;
for(int i=0;i<n;i++) {
scanf("%lf%lf",&list[i].x,&list[i].y);
}
Graham(list,p,n);
n=top;
for(int i=0;i<m;i++) {
scanf("%lf%lf",&list[i].x,&list[i].y);
}
Graham(list,q,m);
m=top;
p[n]=p[0];
q[m]=q[0];
printf("%.5f\n",min(rotating(p,n,q,m),rotating(q,m,p,n)));
}
return 0;
}
/*半平面相交（直线切割多边形）（点标号从1开始）
n^2算法
*/
Point points[MAXN],p[MAXN],q[MAXN];
int n;
double r;
int cCnt,curCnt;
inline void getline(Point x,Point y,double &a,double &b,double &c){
a = y.y - x.y;
b = x.x - y.x;
c = y.x * x.y - x.x * y.y;
}
inline void initial(){
for(int i = 1; i <= n; ++i)p[i] = points[i];
p[n+1] = p[1];
p[0] = p[n];
cCnt = n;
}
inline Point intersect(Point x,Point y,double a,double b,double c){
double u = fabs(a * x.x + b * x.y + c);
double v = fabs(a * y.x + b * y.y + c);
return Point( (x.x * v + y.x * u) / (u + v) , (x.y * v + y.y * u) / (u + v) );
}
inline void cut(double a,double b ,double c){
curCnt = 0;
for(int i = 1; i <= cCnt; ++i){
if(a*p[i].x + b*p[i].y + c >= EPS)q[++curCnt] = p[i];
else {
if(a*p[i-1].x + b*p[i-1].y + c > EPS){
q[++curCnt] = intersect(p[i],p[i-1],a,b,c);
}
if(a*p[i+1].x + b*p[i+1].y + c > EPS){
q[++curCnt] = intersect(p[i],p[i+1],a,b,c);
}
}
}
for(int i = 1; i <= curCnt; ++i)p[i] = q[i];
p[curCnt+1] = q[1];p[0] = p[curCnt];
cCnt = curCnt;
}
inline void solve(){
//注意：默认点是顺时针，如果题目不是顺时针，规整化方向
initial();
for(int i = 1; i <= n; ++i){
double a,b,c;
getline(points[i],points[i+1],a,b,c);
cut(a,b,c);
}
/*
如果要向内推进r，用该部分代替上个函数
for(int i = 1; i <= n; ++i){
Point ta, tb, tt;
tt.x = points[i+1].y - points[i].y;
tt.y = points[i].x - points[i+1].x;
double k = r / sqrt(tt.x * tt.x + tt.y * tt.y);
tt.x = tt.x * k;
tt.y = tt.y * k;
ta.x = points[i].x + tt.x;
ta.y = points[i].y + tt.y;
tb.x = points[i+1].x + tt.x;
tb.y = points[i+1].y + tt.y;
double a,b,c;
getline(ta,tb,a,b,c);
cut(a,b,c);
}
*/
//多边形核的面积
double area = 0;
for(int i = 1; i <= curCnt; ++i)
area += p[i].x * p[i + 1].y - p[i + 1].x * p[i].y;
area = fabs(area / 2.0);
//此时cCnt为最终切割得到的多边形的顶点数，p为存放顶点的数组
}
inline void GuiZhengHua(){
//规整化方向，逆时针变顺时针，顺时针变逆时针
for(int i = 1; i < (n+1)/2; i ++)
swap(points[i], points[n-i]);//头文件加iostream
}
inline void init(){
for(int i = 1; i <= n; ++i)points[i].input();
points[n+1] = points[1];
}
/*
半平面交nlogn算法
*/
#include<stdio.h>
#include<math.h>
#include<algorithm>
#define zero(a) fabs(a)<eps
using namespace std;
const double eps=1e-10;
struct Point {
double x,y;
Point () {};
Point (double xx,double yy) {x=xx;y=yy;};
double operator^(const Point b)const{
return x*b.y-y*b.x;
}
Point operator-(const Point b)const{
return Point(x-b.x,y-b.y);
}
};
struct Line {
Point s,e;
double angle;
Line(){};
Line(Point ss,Point ee){s=ss;e=ee;}
Point operator &(const Line b)const {
Point res=s;
double t=((s-b.e)^(b.s-b.e))/((s-e)^(b.s-b.e));
res.x+=(e.x-s.x)*t;
res.y+=(e.y-s.y)*t;
return res;
}
void puts(){
printf("%.2f %.2f %.2f %.2f %f\n",s.x,s.y,e.x,e.y,angle);
}
};
int n;
Line dq[20005],l[20005];
Point p[20005];
int cmp(Line a,Line b){
if (fabs(a.angle-b.angle)<eps) return ((a.e-a.s)^(b.s-a.s))<-eps;
else return a.angle>b.angle;
}
double xmul(Point p0,Point p1,Point p2){
return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);
}
void HPI(){
sort(l,l+n,cmp);
int ln=0;
for(int i=1,j=0;i<n;i++) {
if (!zero(l[i].angle-l[j].angle)) l[++j]=l[i];
ln=j+1;
}
n=ln;
int bot=0,top=1;
dq[0]=l[0]; dq[1]=l[1];
for(int i=2;i<n;i++){
while(bot<top&&xmul(l[i].s,l[i].e,(dq[top]&dq[top-1]))<-eps) top--;
while(bot<top&&xmul(l[i].s,l[i].e, (dq[bot]&dq[bot+1]))<-eps) bot++;
dq[++top]=l[i];
}
while(bot<top&&xmul(dq[bot].s,dq[bot].e, (dq[top]&dq[top-1]))<-eps) top--;
while(bot<top&&xmul(dq[top].s,dq[top].e, (dq[bot]&dq[bot+1]))<-eps) bot++;
dq[top+1]=dq[bot];
int cnt=0;
for(int i=bot;i<=top;i++) p[cnt++]=dq[i]&dq[i+1];
p[cnt]=p[0];
if (top<=bot+1) cnt=0;
double s=0;
for(int i=0;i<cnt;i++) s+=p[i]^p[i+1];
s=fabs(s*0.5);
printf("%.1f\n",s);
}
int main() {
while(scanf("%d",&n)!=EOF) {
l[0]=Line(Point(0,0),Point(10000,0)); l[n].angle=atan2(l[n].e.y-l[n].s.y,l[n].e.x-l[n].s.x);
l[0+1]=Line(Point(10000,0),Point(10000,10000)); l[0+1].angle=atan2(l[0+1].e.y-l[0+1].s.y,l[0+1].e.x-l[0+1].s.x);
l[0+2]=Line(Point(10000,10000),Point(0,10000)); l[0+2].angle=atan2(l[0+2].e.y-l[0+2].s.y,l[0+2].e.x-l[0+2].s.x);
l[0+3]=Line(Point(0,10000),Point(0,0)); l[0+3].angle=atan2(l[0+3].e.y-l[0+3].s.y,l[0+3].e.x-l[0+3].s.x);
for(int i=0;i<n;i++) {
scanf("%lf%lf%lf%lf",&l[i+4].s.x,&l[i+4].s.y,&l[i+4].e.x,&l[i+4].e.y);
l[i+4].angle=atan2(l[i+4].e.y-l[i+4].s.y,l[i+4].e.x-l[i+4].s.x);
}
n=n+4;
HPI();
}
}
//过三点求圆心坐标
Point waixin(Point a,Point b,Point c){
double a1 = b.x - a.x, b1 = b.y - a.y, c1 = (a1*a1 + b1*b1)/2;
double a2 = c.x - a.x, b2 = c.y - a.y, c2 = (a2*a2 + b2*b2)/2;
double d = a1*b2 - a2*b1;
return Point(a.x + (c1*b2 - c2*b1)/d, a.y + (a1*c2 -a2*c1)/d);
}
//两个圆的公共部分面积
double Area_of_overlap(Point c1,double r1,Point c2,double r2){
double d = dist(c1,c2);
if(r1 + r2 < d + eps)return 0;
if(d < fabs(r1 - r2) + eps){
double r = min(r1,r2);
return PI*r*r;
}
double x = (d*d + r1*r1 - r2*r2)/(2*d);
double t1 = acos(x / r1);
double t2 = acos((d - x)/r2);
return r1*r1*t1 + r2*r2*t2 - d*r1*sin(t1);
}
/*
Pick公式
顶点坐标均是整点的简单多边形：面积=内部格点数目+边上格点数目/2-1
*/
/*
圆的并
*/
/*
题目：CIRU2
题目来源：SPOJ 8119 https://www.spoj.pl/problems/CIRUT/
题目难度：中等偏难
题目内容或思路：
圆的面积并
题意：给n个圆（最多1000个），分别求出覆盖1层的面积、覆盖2层的面积、
覆盖3层的面积。。。覆盖n层的面积
方法见AC大牛blog：
http://hi.baidu.com/aekdycoin/blog/item/c1b28e3711246b3f0b55a95e.html
做题日期：2011.3.27
*/
#include <cstdio>
#include <cstdlib>
#include <climits>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <string>
#include <queue>
#include <map>
#include <vector>
#include <bitset>
#include <cmath>
#include <set>
#include <utility>
#include <ctime>
#define sqr(x) ((x)*(x))
using namespace std;
const int N = 1010;
const double eps = 1e-8;
const double pi = acos(-1.0);
double area[N];
int n;
int dcmp(double x) {
if (x < -eps) return -1; else return x > eps;
}
struct cp {
double x, y, r, angle;
int d;
cp(){}
cp(double xx, double yy, double ang = 0, int t = 0) {
x = xx; y = yy; angle = ang; d = t;
}
void get() {
scanf("%lf%lf%lf", &x, &y, &r);
d = 1;
}
}cir[N], tp[N * 2];
double dis(cp a, cp b) {
return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y));
}
double cross(cp p0, cp p1, cp p2) {
return (p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x);
}
int CirCrossCir(cp p1, double r1, cp p2, double r2, cp &cp1, cp &cp2) {//两圆求交点数 ,并用cp1和cp2存储其交点
double mx = p2.x - p1.x, sx = p2.x + p1.x, mx2 = mx * mx;
double my = p2.y - p1.y, sy = p2.y + p1.y, my2 = my * my;
double sq = mx2 + my2, d = -(sq - sqr(r1 - r2)) * (sq - sqr(r1 + r2));//sq为两圆距离的平方和
if (d + eps < 0) return 0; if (d < eps) d = 0; else d = sqrt(d);
double x = mx * ((r1 + r2) * (r1 - r2) + mx * sx) + sx * my2;
double y = my * ((r1 + r2) * (r1 - r2) + my * sy) + sy * mx2;
double dx = mx * d, dy = my * d; sq *= 2;
cp1.x = (x - dy) / sq; cp1.y = (y + dx) / sq;
cp2.x = (x + dy) / sq; cp2.y = (y - dx) / sq;
if (d > eps) return 2; else return 1;
}
bool circmp(const cp& u, const cp& v) {
return dcmp(u.r - v.r) < 0;
}
bool cmp(const cp& u, const cp& v) {
if (dcmp(u.angle - v.angle)) return u.angle < v.angle;
return u.d > v.d;
}
double calc(cp cir, cp cp1, cp cp2) {
double ans = (cp2.angle - cp1.angle) * sqr(cir.r)
- cross(cir, cp1, cp2) + cross(cp(0, 0), cp1, cp2);
return ans / 2;
}
void CirUnion(cp cir[], int n) {
cp cp1, cp2;
sort(cir, cir + n, circmp);
for (int i = 0; i < n; ++i)
for (int j = i + 1; j < n; ++j)
if (dcmp(dis(cir[i], cir[j]) + cir[i].r - cir[j].r) <= 0)
cir[i].d++;//某圆被覆盖的次数
for (int i = 0; i < n; ++i) {
int tn = 0, cnt = 0;
for (int j = 0; j < n; ++j) {
if (i == j) continue;
if (CirCrossCir(cir[i], cir[i].r, cir[j], cir[j].r,
cp2, cp1) < 2) continue;
cp1.angle = atan2(cp1.y - cir[i].y, cp1.x - cir[i].x);//求圆心指向交点的向量的极角，注意这里atan2(y,x)函数要先y后x
cp2.angle = atan2(cp2.y - cir[i].y, cp2.x - cir[i].x);
//
cp1.d = 1; tp[tn++] = cp1;
cp2.d = -1; tp[tn++] = cp2;
if (dcmp(cp1.angle - cp2.angle) > 0) cnt++;
}
tp[tn++] = cp(cir[i].x - cir[i].r, cir[i].y, pi, -cnt);//第i个圆的左端点
tp[tn++] = cp(cir[i].x - cir[i].r, cir[i].y, -pi, cnt);
sort(tp, tp + tn, cmp);//把交点按极角排序
int p, s = cir[i].d + tp[0].d;
for (int j = 1; j < tn; ++j) {
p = s; s += tp[j].d;
printf("%f %f %f %f\n",tp[j-1].x,tp[j].y,tp[j-1].x,tp[j].y);
area[p] += calc(cir[i], tp[j - 1], tp[j]);//所有被覆盖了p次的圆弧的弦连起来形成的面积，等于覆盖了p次及p次以上的面积，可画图观察
}
}
}
void solve() {
for (int i = 0; i < n; ++i)
cir[i].get();
memset(area, 0, sizeof(area));
CirUnion(cir, n);
for (int i = 1; i <= n; ++i) {
area[i] -= area[i + 1];
printf("[%d] = %.3lf\n", i, area[i]);
}
}
int main() {
while (scanf("%d", &n) != EOF) {
solve();
}
return 0;
}