多边形相交面积模板

  1 /*
  2     类型：多边形相交面积模板
  3 */
  4
  5 #include<cstdio>
  6 #include<iostream>
  7 #include<algorithm>
  8 #include<cstring>
  9 #include<cmath>
 10 using namespace std;
 11 #define maxn 510
 12 const double eps=1E-8;
 13 int sig(double d){
 14     return(d>eps)-(d<-eps);
 15 }
 16 struct Point{
 17     double x,y; Point(){}
 18     Point(double x,double y):x(x),y(y){}
 19     bool operator==(const Point&p)const{
 20         return sig(x-p.x)==0&&sig(y-p.y)==0;
 21     }
 22 };
 23 double cross(Point o,Point a,Point b){
 24     return(a.x-o.x)*(b.y-o.y)-(b.x-o.x)*(a.y-o.y);
 25 }
 26 double area(Point* ps,int n){
 27     ps[n]=ps[0];
 28     double res=0;
 29     for(int i=0;i<n;i++){
 30         res+=ps[i].x*ps[i+1].y-ps[i].y*ps[i+1].x;
 31     }
 32     return res/2.0;
 33 }
 34 int lineCross(Point a,Point b,Point c,Point d,Point&p){
 35     double s1,s2;
 36     s1=cross(a,b,c);
 37     s2=cross(a,b,d);
 38     if(sig(s1)==0&&sig(s2)==0) return 2;
 39     if(sig(s2-s1)==0) return 0;
 40     p.x=(c.x*s2-d.x*s1)/(s2-s1);
 41     p.y=(c.y*s2-d.y*s1)/(s2-s1);
 42     return 1;
 43 }
 44 //多边形切割
 45 //用直线ab切割多边形p，切割后的在向量(a,b)的左侧，并原地保存切割结果
 46 //如果退化为一个点，也会返回去,此时n为1
 47 void polygon_cut(Point*p,int&n,Point a,Point b){
 48     static Point pp[maxn];
 49     int m=0;p[n]=p[0];
 50     for(int i=0;i<n;i++){
 51         if(sig(cross(a,b,p[i]))>0) pp[m++]=p[i];
 52         if(sig(cross(a,b,p[i]))!=sig(cross(a,b,p[i+1])))
 53             lineCross(a,b,p[i],p[i+1],pp[m++]);
 54     }
 55     n=0;
 56     for(int i=0;i<m;i++)
 57         if(!i||!(pp[i]==pp[i-1]))
 58             p[n++]=pp[i];
 59     while(n>1&&p[n-1]==p[0])n--;
 60 }
 61 //---------------华丽的分隔线-----------------//
 62 //返回三角形oab和三角形ocd的有向交面积,o是原点//
 63 double intersectArea(Point a,Point b,Point c,Point d){
 64     Point o(0,0);
 65     int s1=sig(cross(o,a,b));
 66     int s2=sig(cross(o,c,d));
 67     if(s1==0||s2==0)return 0.0;//退化，面积为0
 68     if(s1==-1) swap(a,b);
 69     if(s2==-1) swap(c,d);
 70     Point p[10]={o,a,b};
 71     int n=3;
 72     polygon_cut(p,n,o,c);
 73     polygon_cut(p,n,c,d);
 74     polygon_cut(p,n,d,o);
 75     double res=fabs(area(p,n));
 76     if(s1*s2==-1) res=-res;return res;
 77 }
 78 //求两多边形的交面积
 79 double intersectArea(Point*ps1,int n1,Point*ps2,int n2){
 80     if(area(ps1,n1)<0) reverse(ps1,ps1+n1);
 81     if(area(ps2,n2)<0) reverse(ps2,ps2+n2);
 82     ps1[n1]=ps1[0];
 83     ps2[n2]=ps2[0];
 84     double res=0;
 85     for(int i=0;i<n1;i++){
 86         for(int j=0;j<n2;j++){
 87             res+=intersectArea(ps1[i],ps1[i+1],ps2[j],ps2[j+1]);
 88         }
 89     }
 90     return res;//assumeresispositive!
 91 }
 92 //hdu-3060求两个任意简单多边形的并面积
 93 Point ps1[maxn],ps2[maxn];
 94 int n1,n2;
 95 int main(){
 96     while(scanf("%d%d",&n1,&n2)!=EOF){
 97         for(int i=0;i<n1;i++)
 98             scanf("%lf%lf",&ps1[i].x,&ps1[i].y);
 99         for(int i=0;i<n2;i++)
100             scanf("%lf%lf",&ps2[i].x,&ps2[i].y);
101         double ans=intersectArea(ps1,n1,ps2,n2);
102         ans=fabs(area(ps1,n1))+fabs(area(ps2,n2))-ans;//容斥
103         printf("%.2f\n",ans);
104     }
105     return 0;
106 }