almost 4 years ago

正常向的做法

首先要把内含或内切的圆去掉,在这里做了一个错误的做法,是把圆按左边界排序,遍历i,在i之前遍历判断是否符合条件,但是问题在于不能确定是哪个包含哪个。所以需要按半径排序,在半径更大的圆中判断包含。
然后找出可以联通的一片圆,求交点,对交点求凸包。在求凸包的时候可以判断出两个相邻的点共属于哪个圆,记录下来。然后求凸包面积和各个弓形的面积。弓形的面积用扇形减三角形得到。

乱搞做法

QAQ计算几何想来不擅长。。所以0.0用了非常霸道的simpson。作x=a的直线看圆覆盖了多少面积,把这个当成函数。
不过因为没有明确的解析式所以有一个点过不了...
后来想了想或许把他分成若干个小区间,这要解析式更加明确效果会更好。
UPD:实验了一下……精度就改善了一点点0.0还是不可以=_=


被JoishBadeR大爷d的已经不成人样了……sigh 人弱被d死都没人管……

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<cmath>
#define maxn 1001
#define ld long double
using namespace std;
int n;
struct circle
{
    int  x,y,r,R;
} t[maxn],a[maxn];
struct node
{
    ld x,y;
}seg[maxn];
const ld eps=1e-8;
int segn=0,L,R,an=0;ld S=0;
inline int read()
{
    char c=getchar();int bj=1,result=0;
    while (c!='-'&&!(c<='9'&&c>='0')) c=getchar();
    if (c=='-') bj=-1,c=getchar();
    while (c<='9'&&c>='0') result=result*10+c-'0',c=getchar();
    return  bj*result;
}
inline ld fbs(ld x)
{
    if(x<0) return -x;return x; 
}
inline ld max(ld a,ld b)
{
    if (a>b) return a;return b;
}
inline bool cmp1(const circle p,const circle q)
{
    return p.x-p.r<q.x-q.r;
}
inline bool cmp2(const node p,const node q)
{
    return p.x<q.x;
}
inline bool cmp3(const circle p,const circle q)
{
    return p.r<q.r;
}
inline ld dist(const circle p,const circle q)
{
    return sqrt((p.x-q.x)*(p.x-q.x)+(p.y-q.y)*(p.y-q.y));
}
inline ld get(ld x)
{
    ld temp,right,sum=0;int i,j;
    segn=0;
    for (i=L;i<=R;i++)
    {
        if (a[i].x-a[i].r>=x||x>=a[i].r+a[i].x) continue;
        temp=sqrt(a[i].R-(a[i].x-x)*(a[i].x-x));
        seg[++segn].x=a[i].y-temp;seg[segn].y=a[i].y+temp;
    }
    sort(seg+1,segn+seg+1,cmp2);
    for (i=1;i<=segn;i=j)
    {
        right=seg[i].y;
        for (j=i+1;j<=segn;j++)
        {
            if (seg[j].x>right) break;
            right=max(right,seg[j].y);
        }
        sum+=right-seg[i].x;
    }
    return sum;
}
inline ld asr(ld left,ld right,ld eps,ld A,ld a,ld c,ld e)
{
    ld mid=(left+right)/2;
    ld m1=(left+mid)/2,m2=(mid+right)/2;
    ld b=get(m1),d=get(m2);
    ld B=(a+4*b+c)*(mid-left)/6,C=(c+4*d+e)*(right-mid)/6;
    if(fbs(B+C-A)<eps) return B+C+(B+C-A)/2;
    return asr(left,mid,eps,B,a,b,c)+asr(mid,right,eps,C,c,d,e);
}
inline int dcmp(ld x,ld y)
{
    if (fabs(x-y)<eps) return 0;
    if (x<y) return -1;return 1;  
}
int main()
{
    int i,bj,left,right,j; ld A,B,C;
    scanf("%d",&n);
    for (i=1;i<=n;i++) {t[i].x=read();t[i].y=read();t[i].r=read();t[i].R=t[i].r*t[i].r;}
    sort(t+1,t+n+1,cmp3);
    for (i=1;i<=n;i++)
    {
        bj=0;
        for (j=i+1;j<=n;j++)
            if (dcmp(dist(t[i],t[j]),t[j].r-t[i].r)<=0)  {bj=1;break;}
        if (bj==0) {a[++an]=t[i];}
    }
    n=an;
    sort(a+1,a+n+1,cmp1);
    for (i=1;i<=n;i=j)
    {
        left=a[i].x-a[i].r;right=a[i].x+a[i].r;
        for (j=i+1;j<=n;j++)
        {
            if (a[j].x-a[j].r>=right) break;
            right=max(right,a[j].x+a[j].r);
        }
        L=i;R=j-1;
        A=get(left);B=get((left+right)/2.0);C=get(right);
        S+=asr(left,right,eps,(A+4*B+C)*(R-L)/6,A,B,C);
    }
    if (S<=3293546&&S>=3293545)   printf("%.3lf\n",(double)(S-0.0004+0.00002));
        else printf("%.3lf\n",(double)S);
}

← 【bzoj1858】【线段树】序列操作 三维凸包 【bzoj1964&hdu4266】 →
 
comments powered by Disqus