HDU4130 Shadow

除草,把最近的10k代码拉出来表演

#include<bits/stdc++.h>
#define deb(x)  cerr<< #x <<" = " <<x<<endl
using namespace std;
const int maxn=1010;
const double eps=5e-9;
const double pi=acos(-1);
double form(double x){
    while(x>=2*pi)x-=2*pi;
    while(x<0)x+=2*pi;
    return x;
}
int sgn(double x){return (x>eps)-(x<-eps);}
int sgn(double a,double b){return sgn(a-b);}
double sqr(double x){return x*x;}
namespace Union{

    struct P{
        double x,y;
        P(){}
        P(double x,double y):x(x),y(y){}
        double len2(){
            return sqr(x)+sqr(y);
        }
        double len(){
            return sqrt(len2());
        }
        void print(){
            fprintf(stderr,"(%.3f,%.3f)\n",x,y);
        }
        P turn90(){return P(-y,x);}
        P norm(){return P(x/len(),y/len());}
        P rot(double rad){
            return P(x*cos(rad)-y*sin(rad),x*sin(rad)+y*cos(rad));
        }
        void read(){
            scanf("%lf%lf",&x,&y);
        }
    };
    bool operator<(P a,P b){
        return sgn(a.x,b.x) ? a.x<b.x : a.y<b.y;
    }
    bool operator==(P a,P b){
        return !sgn(a.x-b.x) and !sgn(a.y-b.y);
    }
    P operator+(P a,P b){
        return P(a.x+b.x,a.y+b.y);
    }
    P operator-(P a,P b){
        return P(a.x-b.x,a.y-b.y);
    }
    P operator*(P a,double b){
        return P(a.x*b,a.y*b);
    }
    P operator/(P a,double b){
        return P(a.x/b,a.y/b);
    }
    double operator^(P a,P b){
        return a.x*b.x + a.y*b.y;
    }
    double operator*(P a,P b){
        return a.x*b.y - a.y*b.x;
    }
    double det(P a,P b,P c){
        return (b-a)*(c-a);
    }
    double dis(P a,P b){
        return (b-a).len();
    }
    double Area(vector<P>poly){
        double ans=0;
        for(int i=1;i<poly.size();i++)
            ans+=(poly[i]-poly[0])*(poly[(i+1)%poly.size()]-poly[0]);
        return ans/2;
    }
    struct L{
        P a,b;
        L(){}
        L(P a,P b):a(a),b(b){}
        P v(){return b-a;}
    };
    bool onLine(P p,L l){
        return sgn((l.a-p)*(l.b-p))==0;
    }
    bool onSeg(P p,L s){
        return onLine(p,s) and sgn((s.b-s.a)^(p-s.a))>=0 and sgn((s.a-s.b)^(p-s.b))>=0;
    }
    bool parallel(L l1,L l2){
        return sgn(l1.v()*l2.v())==0;
    }
    P intersect(L l1,L l2){
        double s1=det(l1.a,l1.b,l2.a);
        double s2=det(l1.a,l1.b,l2.b);
        return (l2.a*s2-l2.b*s1)/(s2-s1);
    }
    P project(P p,L l){
        return l.a+l.v()*((p-l.a)^l.v())/l.v().len2();
    }
    double dis(P p,L l){
        return fabs((p-l.a)*l.v())/l.v().len();
    }
    struct C{
        P o;
        double r;
        C(){}
        C(P _o,double _r):o(_o),r(_r){}
        void read(){
            o.read();
            scanf("%lf",&r);
        }
        P at(double rad){
            return o+P(r,0).rot(rad);
        }
        double ang(P p){
            double t=atan2(p.y-o.y,p.x-o.x);
            return form(t);
        }
    };
    bool operator==(C a,C b){
        return a.o==b.o && !sgn(a.r-b.r);
    }
    bool inPoly(P p,vector<P>poly,bool strict){
        int cnt=0;
        for(int i=0;i<poly.size();i++){
            P a=poly[i],b=poly[(i+1)%poly.size()];
            if(onLine(p,L(a,b)))
                return !strict;
            int x=sgn(det(a,p,b));
            int y=sgn(a.y-p.y);
            int z=sgn(b.y-p.y);
            cnt+=(x>0&&y<=0&&z>0);
            cnt-=(x<0&&z<=0&&y>0);
        }
        return cnt;
    }
    bool inC(P p,C c,bool strict){
        double d=dis(p,c.o);
        if(strict){
            return sgn(d,c.r)<0;
        }else{
            return sgn(d,c.r)<=0;
        }
    }
    // 求圆与直线的交点
    //turn90() P(-y,x)
    double lfix(double x){return sgn(x)?x:0;}
    bool intersect(C a, L l, P &p1, P &p2) { 
        double x = ((l.a - a.o)^ (l.b - l.a)),
               y = (l.b - l.a).len2(),
               d = x * x - y * ((l.a - a.o).len2() - a.r * a.r);
        if (sgn(d) < 0) return false;
        d = max(d, 0.0);
        P p = l.a - ((l.b - l.a) * (x / y)), delta = (l.b - l.a) * (sqrt(d) / y);
        p1 = p + delta, p2 = p - delta;
        return true;
    }
    // 求圆与圆的交点,注意调用前要先判定重圆
    bool intersect(C a, C b, P &p1, P &p2) { 
        double s1 = (a.o - b.o).len();
        if (sgn(s1 - a.r - b.r) > 0 || sgn(s1 - fabs(a.r - b.r)) < 0) return false;
        double s2 = (a.r * a.r - b.r * b.r) / s1;
        double aa = (s1 + s2) * 0.5, bb = (s1 - s2) * 0.5;
        P o = (b.o - a.o) * (aa / (aa + bb)) + a.o;
        P delta = (b.o - a.o).norm().turn90() * sqrt(lfix(a.r * a.r - aa * aa));
        p1 = o + delta, p2 = o - delta;
        return true;
    }
    // 求点到圆的切点,按关于点的顺时针方向返回两个点
    bool tang(const C &c, const P &p0, P &p1, P &p2) {
        double x = (p0 - c.o).len2(), d = x - c.r * c.r;
        if (d < eps) return false; // 点在圆上认为没有切点
        P p = (p0 - c.o) * (c.r * c.r / x);
        P delta = ((p0 - c.o) * (-c.r * sqrt(d) / x)).turn90();
        p1 = c.o + p + delta;
        p2 = c.o + p - delta;
        return true;
    }
    // 求圆到圆的外共切线,按关于 c1.o 的顺时针方向返回两条线
    vector<L> extan(const C &c1, const C &c2) {
        vector<L> ret;
        if (sgn(c1.r - c2.r) == 0) {
            P dir = c2.o - c1.o;
            dir = (dir * (c1.r / dir.len())).turn90();
            ret.push_back(L(c1.o + dir, c2.o + dir));
            ret.push_back(L(c1.o - dir, c2.o - dir));
        } else {
            P p = (c1.o * -c2.r + c2.o * c1.r) / (c1.r - c2.r);
            P p1, p2, q1, q2;
            if (tang(c1, p, p1, p2) && tang(c2, p, q1, q2)) {
                if (c1.r < c2.r) swap(p1, p2), swap(q1, q2);
                ret.push_back(L(p1, q1));
                ret.push_back(L(p2, q2));
            }
        }
        return ret;
    }
    // 求圆到圆的内共切线,按关于 c1.o 的顺时针方向返回两条线
    vector<L> intan(const C &c1, const C &c2) {
        vector<L> ret;
        P p = (c1.o * c2.r + c2.o * c1.r) / (c1.r + c2.r);
        P p1, p2, q1, q2;
        if (tang(c1, p, p1, p2) && tang(c2, p, q1, q2)) { // 两圆相切认为没有切线
            ret.push_back(L(p1, q1));
            ret.push_back(L(p2, q2));
        }
        return ret;
    }
    C c[maxn];
    vector<P> poly[maxn];
    int n,m;


    vector<P> intersect(vector<P>poly,C c){
        vector<P>ans;
        for(int i=0;i<poly.size();i++){
            L l=L(poly[i],poly[(i+1)%poly.size()]);
            if(!sgn(l.v().len()))
                continue;
            P p1,p2;
            if(intersect(c,l,p1,p2)){
                if(onSeg(p1,l)){
                    ans.push_back(p1);
                    //cerr<<"ins"<<endl;
                    //p1.print();
                }
                if(onSeg(p2,l))
                    ans.push_back(p2);
            }
        }
        return ans;
    }

    double calcCir(C cir){
        vector<double>ang;
        ang.push_back(0);
        ang.push_back(pi);
        double ans=0;
        for(int i=1;i<=n;i++){
            if(cir==c[i])continue;
            P p1,p2;
            if(intersect(cir,c[i],p1,p2)){
                ang.push_back(form(cir.ang(p1)));
                ang.push_back(form(cir.ang(p2)));
            }
        }
        
        for(int i=1;i<=m;i++){
            vector<P>tmp;
            tmp=intersect(poly[i],cir);
            //deb(tmp.size());
            for(int j=0;j<tmp.size();j++){
            //    tmp[j].print();
                ang.push_back(form(cir.ang(tmp[j])));
            }
        }
        sort(ang.begin(),ang.end());
        for(int i=0;i<ang.size();i++){
            double t1=ang[i],t2=(i+1==ang.size()?ang[0]+2*pi:ang[i+1]);
            P p=cir.at((t1+t2)/2);
            int ok=1;
            for(int j=1;j<=n;j++){
                if(cir==c[j])continue;
                if(inC(p,c[j],true)){
                    ok=0;
                    break;
                }
            }
            for(int j=1;j<=m&&ok;j++){
                if(inPoly(p,poly[j],true)){
                    ok=0;
                    break;
                }
            }
            if(ok){
                double r=cir.r,x0=cir.o.x,y0=cir.o.y;
                ans+=(r*r*(t2-t1)+r*x0*(sin(t2)-sin(t1))-r*y0*(cos(t2)-cos(t1)))/2;        
                /*cerr<<"---"<<endl;
                cerr<<t1<<" "<<t2<<endl;
                cerr<<x0<<" "<<y0<<endl;
                cerr<<ans<<endl;*/
            }
        }
        return ans;
    }
    P st;
    bool bySt(P a,P b){
        return dis(a,st)<dis(b,st);
    }
    double calcSeg(L l){
        double ans=0;
        vector<P>pt;
        pt.push_back(l.a);
        pt.push_back(l.b);

        /*cerr<<"line "<<endl;
        l.a.print();
        l.b.print();
*/
        for(int i=1;i<=n;i++){
            P p1,p2;
            if(intersect(c[i],l,p1,p2)){
                if(onSeg(p1,l))
                    pt.push_back(p1);
                if(onSeg(p2,l))
                    pt.push_back(p2);
            }
        }
        st=l.a;
        sort(pt.begin(),pt.end(),bySt);
        for(int i=0;i+1<pt.size();i++){
            P p1=pt[i],p2=pt[i+1];
            P p=(p1+p2)/2;
            int ok=1;
            for(int j=1;j<=n;j++){
                if(sgn(dis(p,c[j].o),c[j].r)<0){
                    ok=0;
                    break;
                }
            }
            if(ok){
                double x1=p1.x,y1=p1.y,x2=p2.x,y2=p2.y;    
        //        fprintf(stderr,"(%.4f,%.4f)->(%.4f,%.4f)\n",x1,y1,x2,y2);
                double res=(x1*y2-x2*y1)/2;
                //deb(res);
                ans+=res;
            }
        }
        //cerr<<endl;
        return ans;
    }

    C cc[maxn];
    double calc(){

        int nn=0;
        m=1;
        for(int i=1;i<=n;i++){
            int ok=1;
            for(int j=i+1;j<=n;j++){
                if(c[i]==c[j]){
                    ok=0;
                }
            }
            if(ok)
                cc[++nn]=c[i];
        }
        n=nn;
        for(int i=1;i<=n;i++)
            c[i]=cc[i];

//        deb(poly[1].size());
        /*deb(Area(poly[1]));
        for(int i=0;i<poly[1].size();i++)
            poly[1][i].print();*/
        if(Area(poly[1])<0){
            reverse(poly[1].begin(),poly[1].end());
        }
        double ans=0;

        for(int i=1;i<=n;i++){
            ans+=calcCir(c[i]);
            //cerr<<sqr(c[i].r)*pi<<endl;
        }
        double res=0;
        for(int i=0;i<poly[1].size();i++){
            P p1=poly[1][i],p2=poly[1][(i+1)%poly[1].size()];
            if(p1==p2)continue;
            res+=calcSeg(L(p1,p2));
        }
    //    cerr<<res<<endl;
//        cerr<<"Area"<<endl;
        ans+=res;
        return ans;
    }
    vector<P> convex(vector<P>p){
        sort(p.begin(),p.end());
        /*deb(p.size());
        for(int i=0;i<p.size();i++)
            p[i].print();*/
        vector<P>ans,S;
        for(int i=0;i<p.size();i++){
            while(S.size()>=2
                    && sgn(det(S[S.size()-2],S.back(),p[i]))<=0)
                S.pop_back();
            S.push_back(p[i]);
        }//dw
        ans=S;
        S.clear();
        for(int i=(int)p.size()-1;i>=0;i--){
            while(S.size()>=2
                    && sgn(det(S[S.size()-2],S.back(),p[i]))<=0)
                S.pop_back();
            S.push_back(p[i]);        
        }//up    
        for(int i=1;i+1<S.size();i++)
            ans.push_back(S[i]);
        /*deb(ans.size());
        for(int i=0;i<ans.size();i++)
            ans[i].print();
        exit(0);*/
        return ans;
    }
}


struct P{
    double x,y,z;
    P(){}
    P(double _x,double _y,double _z):x(_x),y(_y),z(_z){}
    void read(){
        scanf("%lf%lf%lf",&x,&y,&z);
    }
    double len(){
        return sqrt(x*x+y*y+z*z);
    }
}p[maxn];
double r[maxn];
double operator^(P a,P b){
    return a.x*b.x+a.y*b.y+a.z*b.z;
}
P operator*(P a,P b){
    return P(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x);
}
P operator/(P a,double b){
    return P(a.x/b,a.y/b,a.z/b);
}


int main(){
    int n,m;
    P D;
    while(cin>>n>>m){
        Union::n=n;
        D.read();
        if(!n&&!m&&!D.x&&!D.y&&!D.z)break;
        for(int i=1;i<=n;i++){
            p[i].read();
            scanf("%lf",&r[i]);
        }
        D=D/D.len();
        P X,Y,Z=D;
        if(sgn((P(1,0,0)*Z).len())){
            X=P(1,0,0)*Z;
            Y=Z*X;
        }else{    
            X=P(0,1,0)*Z;
            Y=Z*X;
        }
        X=X/X.len();
        Y=Y/Y.len();
        double rad=acos(Z^P(0,0,-1));
        double rat=cos(rad);
        for(int i=1;i<=n;i++){
            Union::c[i].o.x=(X^p[i]);
            Union::c[i].o.y=(Y^p[i]);
            Union::c[i].r=r[i];
            //Union::c[i].o.print();
            //cerr<<Union::c[i].r<<endl;
        }
        //cerr<<endl;
        vector<Union::P>pt,ch;
        pt.resize(m);
        for(int i=0;i<m;i++){
            P p;p.read();
            pt[i].x=(X^p);
            pt[i].y=(Y^p);
        }
        ch=Union::convex(pt);
        //for(int i=0;i<ch.size();i++)
    //        ch[i].print();
        Union::poly[1]=ch;
    //    cerr<<Union::calc()<<endl;
        printf("%.4f\n",Union::calc()/rat);
    }
    return 0;
}

发表回复

您的电子邮箱地址不会被公开。 必填项已用*标注