zoukankan      html  css  js  c++  java
  • 计算几何+求质心+求多边形与圆交面积——ICPC GNYR 2019

    /*
    先求出多边形的质心
        将多边形分割成三角形 OViVi+1,
        求每个三角形的质心((x1+x2+x3)/3,(y1+y2+y3)/3),然后再有向面积加权 
    然后求圆和多边形的交点:套模板即可 
    */
    #include<bits/stdc++.h>
    using namespace std;
    
    typedef double db;
    const db eps=1e-6;
    const db pi=acos(-1);
    int sign(db k){
        if (k>eps) return 1; else if (k<-eps) return -1; return 0;
    }
    int cmp(db k1,db k2){return sign(k1-k2);}
    int inmid(db k1,db k2,db k3){return sign(k1-k3)*sign(k2-k3)<=0;}// k3 在 [k1,k2] 内 
    struct point{
        db x,y;
        point operator + (const point &k1) const{return (point){k1.x+x,k1.y+y};}
        point operator - (const point &k1) const{return (point){x-k1.x,y-k1.y};}
        point operator * (db k1) const{return (point){x*k1,y*k1};}
        point operator / (db k1) const{return (point){x/k1,y/k1};}
        int operator == (const point &k1) const{return cmp(x,k1.x)==0&&cmp(y,k1.y)==0;}
        point turn(db k1){return (point){x*cos(k1)-y*sin(k1),x*sin(k1)+y*cos(k1)};}
        point turn90(){return (point){-y,x};}
        bool operator < (const point k1) const{
            int a=cmp(x,k1.x);
            if (a==-1) return 1; else if (a==1) return 0; else return cmp(y,k1.y)==-1;
        }
        db abs(){return sqrt(x*x+y*y);}
        db abs2(){return x*x+y*y;}
        db dis(point k1){return ((*this)-k1).abs();}
        point unit(){db w=abs(); return (point){x/w,y/w};}
    };
    int inmid(point k1,point k2,point k3){return inmid(k1.x,k2.x,k3.x)&&inmid(k1.y,k2.y,k3.y);}
    db cross(point k1,point k2){return k1.x*k2.y-k1.y*k2.x;}
    db dot(point k1,point k2){return k1.x*k2.x+k1.y*k2.y;}
    db rad(point k1,point k2){return atan2(cross(k1,k2),dot(k1,k2));}
    point proj(point k1,point k2,point q){ // q 到直线 k1,k2 的投影 
        point k=k2-k1; return k1+k*(dot(q-k1,k)/k.abs2());
    }
    db disSP(point k1,point k2,point q){
        point k3=proj(k1,k2,q);
        if (inmid(k1,k2,k3)) return q.dis(k3); else return min(q.dis(k1),q.dis(k2));
    }
    struct circle{
        point o; db r;
        int inside(point k){return cmp(r,o.dis(k));}
    };
    vector<point> getCL(circle k1,point k2,point k3){ // 求直线和圆的交点,沿着 k2->k3 方向给出 , 相切给出两个 
        point k=proj(k2,k3,k1.o); db d=k1.r*k1.r-(k-k1.o).abs2();
        if (sign(d)==-1) return {};
        point del=(k3-k2).unit()*sqrt(max((db)0.0,d)); return {k-del,k+del};
    }
    db getarea(circle k1,point k2,point k3){
        // 圆 k1 与三角形 k2 k3 k1.o 的有向面积交
        point k=k1.o; k1.o=k1.o-k; k2=k2-k; k3=k3-k;
        int pd1=k1.inside(k2),pd2=k1.inside(k3); 
        vector<point>A=getCL(k1,k2,k3);
        if (pd1>=0){
            if (pd2>=0) return cross(k2,k3)/2;
            return k1.r*k1.r*rad(A[1],k3)/2+cross(k2,A[1])/2;
        } else if (pd2>=0){ 
            return k1.r*k1.r*rad(k2,A[0])/2+cross(A[0],k3)/2;
        }else {
            int pd=cmp(k1.r,disSP(k2,k3,k1.o));
            if (pd<=0) return k1.r*k1.r*rad(k2,k3)/2;
            return cross(A[0],A[1])/2+k1.r*k1.r*(rad(k2,A[0])+rad(A[1],k3))/2;
        }
    }
    
    
    point p[200],k;
    circle c;
    db X,Y,sum;
    int n;
    
    int main(){
        cin>>n;
        for(int i=0;i<n;i++)cin>>p[i].x>>p[i].y;
        k.x=k.y=0;
        for(int i=0;i<n;i++)//求总面积 
            sum+=cross(p[i]-k,p[(i+1)%n]-k)/2;
        for(int i=0;i<n;i++){
            db area=cross(p[i]-k,p[(i+1)%n]-k)/2;
            db x=(p[i].x+p[(i+1)%n].x)/3;
            db y=(p[i].y+p[(i+1)%n].y)/3;
            X+=x*area/sum;
            Y+=y*area/sum;
        }
        c.o=(point){X,Y};
        c.r=sqrt(sum/pi); //sum=pi*r*r
        
        db ans=0; 
        for(int i=0;i<n;i++){ //求相交面积 
            ans+=getarea(c,p[i],p[(i+1)%n]);
        }
        
        printf("%.4lf
    ",ans/sum);
    }
  • 相关阅读:
    纪中第三天
    纪中第一天
    图片验证码的实现
    使用监听器解决路径问题
    log4j测试示例
    redis示例
    kafka示例
    CSRF verification failed. Request aborted.
    TemplateDoesNotExist
    创建 django 项目命令
  • 原文地址:https://www.cnblogs.com/zsben991126/p/12791357.html
Copyright © 2011-2022 走看看