算法几何 ChaoYeon

Chae Yeon is a popular pop female singer who rose to fame with her amazing sexy dance style and the sounds of nature voice she has. She proved to be a great dancer, and she showed off her vocals in her live performances. If you had ever seen her dance, I bet you’d love it.
I felt stage lighting interesting when I was enjoying Chae Yeon’s performance. We all know that stage lighting instruments are used for the concerts and other performances taking place in live performance venues. They are also used to light the stages. Actually this is a color mixing process. There are two types of color mixing: Additive and Subtractive. Most stages use the former and in this case there are three primary colors: red, green, and blue. In the absence of color, or when no colors are showing, the stage is black. If all three primary colors are showing, the result is white. When red and green combine together, the result is yellow. When red and blue combine together, the result is magenta. When blue and green combine together, the result is cyan. When two same color combine together, the result is still that color.

We have got the coordinate and the primary color of the figure that each Stage Lighting Instrument sent out. For simplicity’s sake, we can just treat the figure as a circle. Of course we’ll know the radius of each colored circle. Some color may be changed based on the Color Mixed Theory we mentioned above. Can you find the area of each color?

Input

The first line consists of an integer T, indicating the number of test cases.
The first line of each case consists of three integers R, G, B, indicating the number of red circles, green circles and blue circles. The next R + G + B lines, each line consists of three integer x, y, r, indicating the coordinate and the radius. The first R lines descript the red circles, the second G lines descript the green circles and the last B lines descript the blue circles.

Output

Output seven floating numbers, they are the area of red, green, blue, white, yellow, magenta and cyan. Please take each number with two factional digits.
Constraints
0 < T <= 20
0 <= R, G, B <= 100
-100 <= x, y <= 100; 0 <= r <= 100

Sample Input

1

1 1 1

0 2 3

2 0 3

-2 0 3

Sample Output

9.28 15.04 15.04 4.92 7.04 7.04 1.28

Hint

Gentleman’s Reminder: please make sure that your program won’t output “-0.00”.

Source

The 34th ACM/ICPC Asia Regional Wuhan Site —— Online Contest

 

解题思路

根据容斥原理，S圆1,2交=S圆1,2并-S圆1-S圆2。此原理可推广到两类圆的情况，即：S红,绿交=S红,绿并-S红-S绿。
也能推广到三种颜色圆的情况，即：S红,绿,蓝交=S红,绿,蓝并-S红-S绿-S蓝+S红,绿并+S红,蓝并+S绿,蓝并。
所以本题就转化成了求圆的并的问题，需要求：S红、S绿、S蓝、S红,绿并、S红,蓝并、S绿,蓝并、S红,绿,蓝并。
圆的并的解法：
将每一个圆与其他所有圆的交点找到，并按照与圆心构成的夹角大小来排序，之后逐一取相邻的两点构成圆弧，循环判断圆弧中点是否在其他圆内，如果是，不计算，如果未被覆盖，计算弓形面积+三角形面积。
其中，多个三角形面积构成多边形面积，取P(0,0)，计算公式：(x1*y2-y1*x2)/2。
多注意浮点数的处理。

范例的解法是采用C++实现的。

优化方案一、改用C语言实现，优化排序的算法，原为C++的sort算法，首先不记录交点，改为记录与圆心夹角，采用插入排序，新增一个处理一个。
优化效果：NOJ：90ms—>60ms。
优化方案二、取消排序，以及后面的判断是否被覆盖，初始化为0~2pi，当作区间来处理，在计算出夹角后，每次割掉相交的弧度范围，如果跨越了2pi，则按照两段来处理：Ab~2pi和0~Ae。

优化效果：NOJ：60ms—>20ms。

范例程序C++

/**
 * File Name: ChaeYeon.cpp
 * Created Time:  7/4/2009 5:04:19 PM
 * Author: momodi
 * Description: 
 */
#include <cstdio>
#include <cstring>
#include <cmath>
#include <cstdlib>
#include <algorithm>
#include <vector>
using namespace std;
#define NEXT(v, n) (((v) == (n))? 0 : (v))
#define SZ(x) ((int) x.size())
template <class T> bool get_max(T& a, const T &b) {return b > a? a = b, 1: 0;}
template <class T> bool get_min(T& a, const T &b) {return b < a? a = b, 1: 0;}
#define SQR(v) ((v) * (v))
const double eps = 1e-9;
const double pi = acos(-1.0);

int sgn(const double &a) {
    return (a > eps) - (a < -eps);
}

const int zx[] = {
    0, 1, 0, -1
};
const int zy[] = {
    1, 0, -1, 0
};

struct P {
    double x, y;
    P(const double &_x, const double &_y)
        :x(_x), y(_y) {}
    P() {}
    bool operator == (const P &a) const {
        return sgn(x - a.x) == 0 && sgn(y - a.y) == 0;
    }
    P operator + (const P &a) const {
        return P(x + a.x, y + a.y);
    }
    P operator - (const P &a) const {
        return P(x - a.x, y - a.y);
    }
    P operator * (const double &a) const {
        return P(x * a, y * a);
    }
    P operator / (const double &a) const {
        return P(x / a, y / a);
    }
    P trunc(double a) const {
        a /= sqrt(SQR(x) + SQR(y));
        return P(x * a, y * a);
    }
    P turn_left() const {
        return P(-y, x);
    }
    P turn_right() const {
        return P(y, -x);
    }
    const P& input() {
        scanf("%lf %lf", &x, &y);
        return *this;
    }
    const P& output() const {
        printf("P: %.12lf %.12lf\n", x, y);
        return *this;
    }
};

double dist2(const P &a, const P &b) {
    return SQR(a.x - b.x) + SQR(a.y - b.y);
}
double dist(const P &a, const P &b) {
    return sqrt(SQR(a.x - b.x) + SQR(a.y - b.y));
}
double cross(const P &a, const P &b, const P &c) {
    return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
double dmul(const P &a, const P &b, const P &c) {
    return (b.x - a.x) * (c.x - a.x) + (b.y - a.y) * (c.y - a.y);
}

struct C {
    P mid;
    double r;
    C(const P &_mid, const double &_r)
        :mid(_mid), r(_r) {}
    C() {}
    bool operator == (const C &a) const {
        return mid == a.mid && sgn(r - a.r) == 0;
    }
    bool in(const C &a) const {
        return sgn(r + dist(mid, a.mid) - a.r) <= 0;
    }
    const C &input() {
        mid.input();
        scanf("%lf", &r);
        return *this;
    }
    const C &output() const {
        printf("P: %.12lf %.12lf R: %.12lf\n", mid.x, mid.y, r);
    }
};

double cal_angle(const C &c, const P &a, const P &b) {
    double k = dmul(c.mid, a, b) / SQR(c.r);
    get_min(k, 1.0);
    get_max(k, -1.0);
    return acos(k);
}

double cal_area(const C &c, const P &a, const P &b) {
    return SQR(c.r) * cal_angle(c, a, b) / 2 - cross(c.mid, a, b) / 2;
}

struct cmp {
    P mid;
    cmp(const P &_mid)
        :mid(_mid) {}
    bool operator () (const P &a, const P &b) {
        return atan2(a.y - mid.y, a.x - mid.x) < atan2(b.y - mid.y, b.x - mid.x);
    }
};

bool circles_intersection(const C &a, const C &b, P &c1, P &c2) {
    double dd = dist(a.mid, b.mid);
    if (sgn(dd - (a.r + b.r)) >= 0) {
        return false;
    }
    double l = (dd + (SQR(a.r) - SQR(b.r)) / dd) / 2;
    double h = sqrt(SQR(a.r) - SQR(l));
    c1 = a.mid + (b.mid - a.mid).trunc(l) + (b.mid - a.mid).turn_left().trunc(h);
    c2 = a.mid + (b.mid - a.mid).trunc(l) + (b.mid - a.mid).turn_right().trunc(h);
    return true;
}

bool cover(const C &c, const P &a, const P &b, const vector<C> &cir) {
    P p = c.mid + ((a + b) / 2 - c.mid).trunc(c.r);
    for (vector<C>::const_iterator it = cir.begin(); it != cir.end(); ++it) {
        if (sgn(dist2(p, it->mid) - SQR(it->r)) < 0) {
            return true;
        }
    }
    return false;
}

double cal_area(const vector<C> &in) {
    vector<C> cir;
    for (int i = 0; i < SZ(in); ++i) {
        if (sgn(in[i].r) == 0) {
            continue;
        }
        bool flag = false;
        for (int j = i + 1; j < SZ(in); ++j) {
            if (in[i] == in[j]) {
                flag = true;
                break;
            }
        }
        if (flag) {
            continue;
        }
        for (j = 0; j < SZ(in); ++j) {
            if (!(in[i] == in[j]) && in[i].in(in[j])) {
                flag = true;
                break;
            }
        }
        if (flag) {
            continue;
        }
        cir.push_back(in[i]);
    }
    vector<vector<P> > point_on_circle(SZ(cir));
    for (i = 0; i < SZ(cir); ++i) {
        for (int z = 0; z < 4; ++z) {
            point_on_circle[i].push_back(cir[i].mid + P(zx[z], zy[z]).trunc(cir[i].r));
        }
    }
    for (i = 0; i < SZ(cir); ++i) {
        for (int j = i + 1; j < SZ(cir); ++j) {
            P a, b;
            if (circles_intersection(cir[i], cir[j], a, b)) {
                point_on_circle[i].push_back(a);
                point_on_circle[i].push_back(b);
                point_on_circle[j].push_back(a);
                point_on_circle[j].push_back(b);
            }
        }
    }
    for (i = 0; i < SZ(cir); ++i) {
        sort(point_on_circle[i].begin(), point_on_circle[i].end(), cmp(cir[i].mid));
        point_on_circle[i].erase(unique(point_on_circle[i].begin(), point_on_circle[i].end()), point_on_circle[i].end());
    }
    double ans = 0;
    for (i = 0; i < SZ(cir); ++i) {
        for (int j = 0; j < SZ(point_on_circle[i]); ++j) {
            const P &a = point_on_circle[i][j];
            const P &b = point_on_circle[i][NEXT(j + 1, SZ(point_on_circle[i]))];
            if (!cover(cir[i], a, b, cir)) {
                ans += cross(P(0, 0), a, b) / 2;
                ans += cal_area(cir[i], a, b);
            }
        }
    }
    return ans;
}
double mabs(double a) {
    return a + eps;
}

int main() {
    int ca;
    scanf("%d", &ca);
    while (ca--) {
        int R, G, B;
        scanf("%d %d %d", &R, &G, &B);
        vector<C> r, g, b, rg, rb, gb, rgb;
        for (int i = 0; i < R; ++i) {
            C tmp;
            tmp.input();
            r.push_back(tmp);
            rg.push_back(tmp);
            rb.push_back(tmp);
            rgb.push_back(tmp);
        }
        for (i = 0; i < G; ++i) {
            C tmp;
            tmp.input();
            g.push_back(tmp);
            rg.push_back(tmp);
            gb.push_back(tmp);
            rgb.push_back(tmp);
        }
        for (i = 0; i < B; ++i) {
            C tmp;
            tmp.input();
            b.push_back(tmp);
            rb.push_back(tmp);
            gb.push_back(tmp);
            rgb.push_back(tmp);
        }
        double rs = cal_area(r);
        double gs = cal_area(g);
        double bs = cal_area(b);
        double rgs = rs + gs - cal_area(rg);
        double rbs = rs + bs - cal_area(rb);
        double gbs = gs + bs - cal_area(gb);
        double rgbs = cal_area(rgb) - rs - gs - bs + rgs + rbs + gbs;
        printf("%.2f %.2f %.2f %.2f %.2f %.2f %.2f\n", 
                mabs(rs - rgs - rbs + rgbs), 
                mabs(gs - rgs - gbs + rgbs),
                mabs(bs - rbs - gbs + rgbs),
                mabs(rgbs),
                mabs(rgs - rgbs),
                mabs(rbs - rgbs),
                mabs(gbs - rgbs));
    }
    return 0;
}

优化方案一

#include <stdio.h>
#include "string.h"
#include "math.h"

#define S(x) (x)*(x)
const double pi=3.141592653589,e=1e-9;

typedef struct{
	double x,y,r;
}C;

int Ne(const double a)
{
	return (a>e)-(a<-e);
}

double dis(const double a,const double b,const double c,const double d)
{
	return sqrt(S(a-c)+S(b-d));
}

void InQue(double p[],int *t)
{
	int i,j;
	double tp;
	if(p[*t-1]<0)p[*t-1]+=2*pi;
	if(p[*t-1]>2*pi)p[*t-1]-=2*pi;
	for (i=0;i<*t-1;i++)
	{
		if(p[*t-1]<p[i]){
			tp=p[*t-1];
			for (j=*t-1;j>i;j--)
				p[j]=p[j-1];
			p[i]=tp;
		}
		else if(!Ne(p[*t-1]-p[i])){
			(*t)--;break;}
	}
}

double G_ar(C c[],int k)
{
	char cov;
	int i,j,t=0,q;
	C Z[300];
	double p[600],f=0,d,h,a,a1,a2;
	for (i=0;i<k;i++)
	{
		if(c[i].r==0)continue;
		cov=0;
		for (j=i+1;j<k;j++)
		{
			if(c[i].x==c[j].x && c[i].y==c[j].y && c[i].r==c[j].r){
				cov=1;break;}
		}
		if(cov)continue;
		for (j=0;j<k;j++)
		{
			if (i!=j && c[i].r<c[j].r)
				if (dis(c[i].x,c[i].y,c[j].x,c[j].y)<(c[j].r-c[i].r)+e){
					cov=1;break;}
		}
		if(cov)continue;
		Z[t++]=c[i];
	}

	k=t;
	for (i=0;i<k;i++)
	{
		for (t=j=0;j<k;j++)
		{
			if (i!=j){
				d=dis(Z[i].x,Z[i].y,Z[j].x,Z[j].y);
				if(d<(Z[j].r+Z[i].r)){
					a=(S(d)+S(Z[i].r)-S(Z[j].r))/2/d/Z[i].r;
					a>1 &&	(a=1);
					a<-1 && (a=-1);
					a=acos(a);
					h=atan2(Z[j].y-Z[i].y,Z[j].x-Z[i].x);
					p[t++]=h-a;
					InQue(p,&t);
					p[t++]=h+a;
					InQue(p,&t);
				}
			}
		}

		if (t){
			for (j=0;j<t;j++)
			{
				a1=p[j],a2=p[(j==t-1)?0:(j+1)];
				if(a2<a1)a2+=2*pi;
				d=(a1+a2)/2;
				a=Z[i].x+Z[i].r*cos(d),h=Z[i].y+Z[i].r*sin(d);
				for (q=0;q<k;q++)
					if(i!=q && dis(a,h,Z[q].x,Z[q].y)<Z[q].r)break;
				if(q==k){
					f+=S(Z[i].r)*(a2-a1-sin(a2-a1))/2;
					f+=((Z[i].x+Z[i].r*cos(a1))*(Z[i].y+Z[i].r*sin(a2))-
						(Z[i].x+Z[i].r*cos(a2))*(Z[i].y+Z[i].r*sin(a1)))/2;
				}
			}
			
		}
		else{
			f+=pi*S(Z[i].r);}
	}
	return f;
}

int main()
{
	int T,R,G,B,i;
	C c[300],tc[200];
	double f[7];
	scanf("%d",&T);
	for (;T;T--)
	{
		memset(f,0,56);
		scanf("%d %d %d",&R,&G,&B);
		for (i=0;i<R;i++)
			scanf("%lf %lf %lf",&c[i].x,&c[i].y,&c[i].r);
		f[0]=G_ar(c,R);
		for (i=R;i<R+G;i++)
			scanf("%lf %lf %lf",&c[i].x,&c[i].y,&c[i].r);
		f[3]=G_ar(c,R+G);
		for (i=R+G;i<R+G+B;i++)
			scanf("%lf %lf %lf",&c[i].x,&c[i].y,&c[i].r);
		f[6]=G_ar(c,R+G+B);		
		memcpy(tc,c+R,G*24);
		f[1]=G_ar(tc,G);
		memcpy(tc,c+R+G,B*24);
		f[2]=G_ar(tc,B);
		memcpy(tc,c+R,(G+B)*24);
		f[5]=G_ar(tc,G+B);
		memcpy(c+R,c+(R+G),B*24);
		f[4]=G_ar(c,R+B);
		
		f[3]=f[0]+f[1]-f[3];
		f[4]=f[0]+f[2]-f[4];
		f[5]=f[1]+f[2]-f[5];
		f[6]=f[6]-f[0]-f[1]-f[2]+f[3]+f[4]+f[5];
		printf("%.2f %.2f %.2f %.2f %.2f %.2f %.2f\n",f[0]-f[3]-f[4]+f[6]+e,
			f[1]-f[3]-f[5]+f[6]+e,f[2]-f[4]-f[5]+f[6]+e,f[6]+e,f[3]-f[6]+e,f[4]-f[6]+e,f[5]-f[6]+e);
	}
}

优化方案二

#include <stdio.h>
#include "string.h"
#include "math.h"

#define S(x) (x)*(x)
const double pi=3.141592653589,e=1e-9;

typedef struct{
	double x,y,r;
}C;

typedef struct{
	double b,e;
}A;

double dis(const double a,const double b,const double c,const double d){
	return sqrt(S(a-c)+S(b-d));}

void InQue(A p[],int *t,A nt)
{
	int i,j;
	for (i=0;i<*t;i++)
	{
		if(nt.b>p[i].e)continue;
		else if(nt.b>p[i].b){
			if (nt.e<p[i].e){
				for (j=*t;j>i+1;j--)
					p[j]=p[j-1];
				p[i+1].b=nt.e;
				p[i+1].e=p[i].e;
				p[i].e=nt.b;
				(*t)++;}
			else{
				p[i].e=nt.b;}
		}
		else{
			if (nt.e>p[i].e-e){
				for (j=i;j<*t-1;j++)
					p[j]=p[j+1];
				(*t)--,i--;}
			else if(nt.e>p[i].b){
				p[i].b=nt.e;}
		}
	}
}

double G_ar(C c[],int k)
{
	char cov;
	int i,j,t=0;
	C Z[300];
	A p[10],nt;
	double f=0,d,h,a;
	for (i=0;i<k;i++)
	{
		if(c[i].r==0)continue;
		cov=0;
		for (j=i+1;j<k;j++)
		{
			if(c[i].x==c[j].x && c[i].y==c[j].y && c[i].r==c[j].r){
				cov=1;break;}
		}
		if(cov)continue;
		for (j=0;j<k;j++)
		{
			if (i!=j && c[i].r<c[j].r)
				if (dis(c[i].x,c[i].y,c[j].x,c[j].y)<(c[j].r-c[i].r)+e){
					cov=1;break;}
		}
		if(cov)continue;
		Z[t++]=c[i];
	}

	k=t;
	for (i=0;i<k;i++)
	{
		t=1,p[0].b=0,p[0].e=2*pi;
		for (j=0;j<k;j++)
		{
			if (i!=j){
				d=dis(Z[i].x,Z[i].y,Z[j].x,Z[j].y);
				if(d<(Z[j].r+Z[i].r)){
					a=(S(d)+S(Z[i].r)-S(Z[j].r))/2/d/Z[i].r;
					a>1 && (a=1);
					a<-1 && (a=-1);
					a=acos(a);
					h=atan2(Z[j].y-Z[i].y,Z[j].x-Z[i].x);
					if (h-a<-e && h+a>e)
					{
						nt.b=h-a+2*pi,nt.e=2*pi;
						InQue(p,&t,nt);
						nt.b=0,nt.e=h+a;
						InQue(p,&t,nt);
					}
					else
					{
						nt.b=h-a,nt.b<-e && (nt.b+=2*pi);
						nt.e=h+a,nt.e<e && (nt.e+=2*pi);
						InQue(p,&t,nt);
					}
				}
			}
		}

		for (j=0;j<t;j++)
		{
			a=p[j].b,h=p[j].e;
			f+=S(Z[i].r)*(h-a-sin(h-a))/2;
			f+=((Z[i].x+Z[i].r*cos(a))*(Z[i].y+Z[i].r*sin(h))-
				(Z[i].x+Z[i].r*cos(h))*(Z[i].y+Z[i].r*sin(a)))/2;
		}
	}
	return f;
}

int main()
{
	int T,R,G,B,i;
	C c[300];
	double f[7];
	scanf("%d",&T);
	for (;T;T--)
	{
		memset(f,0,56);
		scanf("%d %d %d",&R,&G,&B);
		for (i=0;i<R;i++)
			scanf("%lf %lf %lf",&c[i].x,&c[i].y,&c[i].r);
		f[0]=G_ar(c,R);
		for (i=R;i<R+G;i++)
			scanf("%lf %lf %lf",&c[i].x,&c[i].y,&c[i].r);
		f[3]=G_ar(c,R+G);
		for (i=R+G;i<R+G+B;i++)
			scanf("%lf %lf %lf",&c[i].x,&c[i].y,&c[i].r);
		f[6]=G_ar(c,R+G+B);
		f[1]=G_ar(c+R,G);
		f[2]=G_ar(c+R+G,B);
		f[5]=G_ar(c+R,G+B);
		memcpy(c+R,c+(R+G),B*24);
		f[4]=G_ar(c,R+B);
		
		f[3]=f[0]+f[1]-f[3];
		f[4]=f[0]+f[2]-f[4];
		f[5]=f[1]+f[2]-f[5];
		f[6]=f[6]-f[0]-f[1]-f[2]+f[3]+f[4]+f[5];
		printf("%.2f %.2f %.2f %.2f %.2f %.2f %.2f\n",f[0]-f[3]-f[4]+f[6]+e,
			f[1]-f[3]-f[5]+f[6]+e,f[2]-f[4]-f[5]+f[6]+e,f[6]+e,f[3]-f[6]+e,f[4]-f[6]+e,f[5]-f[6]+e);
	}
}

调试经验

优化方案二在f[6]崩溃
崩溃原因在浮点数等于的处理不对，下面这句话
     else
     {
      nt.b=h-a,nt.b<-e && (nt.b+=2*pi);
      nt.e=h+a,nt.e<-e && (nt.e+=2*pi);
      InQue(p,&t,nt);
     }
改为下面，表示
     else
     {
      nt.b=h-a,nt.b<-e && (nt.b+=2*pi);
      nt.e=h+a,nt.e<e && (nt.e+=2*pi);
      InQue(p,&t,nt);
     }