半平面交,注意直线方向!!!
对于凸包上任意一条边$LINE(p_i,p_{i+1})$都有$S_{Delta{p_i} {p_{i + 1}}p} < S_{Delta{p_0} {p_1}p}$
如果我们用叉积来算面积:
$P=(x,y)$
$A=p_0=(x_1,y_1)$
$B=p_1=(x_2,y_2)$
$C=p_{i+1}=(x_3,y_3)$(至于为什么是C为i+1而不是D为i+1,画画图就知道了)
$D=p_i=(x_4,y_4)$
就有不等式:
$(x-x_2,y-y_2)(x-x_1,y-y_1)leqslant(x-x_3,y-y_3)(x-x_4,y-y_4)$
$Rightarrow$
$(x-x_2)(y-y_1)-(y-y_2)(x-x_1) leqslant (x-x_3)(y-y_4)-(y-y_3)(x-x_4)$
$Rightarrow$
$(-y_1+y_2-y_3-y_4)x+(x_1-x_2+x_3-x_4)y+$$(x_2 imes y_1-y_2 imes x_1)+(x_4 imes y_3-y_4 imes x_3) leqslant 0$
$Rightarrow$
$(-y_1+y_2-y_3-y_4)x+(x_1-x_2+x_3-x_4)y+overrightarrow B imesoverrightarrow A+overrightarrow D imes overrightarrow Cleqslant 0$
就得到了一个半平面
交一交即可
1 #include<iostream> 2 #include<cstdio> 3 #include<cstdlib> 4 #include<algorithm> 5 #include<cmath> 6 #include<string> 7 #include<cstring> 8 #include<set> 9 #include<queue> 10 #define op operator 11 #define re(i,l,r) for(int i=(l);i<=(r);i++) 12 #define rre(i,r,l) for(int i=(r);i>=(l);i--) 13 using namespace std; 14 template <typename Q> 15 void inin(Q &ret) 16 { 17 ret=0;int f=0;char ch=getchar(); 18 while(ch<'0'||ch>'9'){if(ch=='-')f=1;ch=getchar();} 19 while(ch>='0'&&ch<='9')ret=(ret<<3)+(ret<<1)+ch-'0',ch=getchar(); 20 ret=f?-ret:ret; 21 } 22 int shurux,shuruy; 23 const double eps=1e-8; 24 int dcmp(const double &x){return abs(x)<eps?0:x<0?-1:1;} 25 const double pi=acos(-1); 26 const double f(const double &a){return a*a;} 27 struct xl 28 { 29 double x,y; 30 xl(double x=0,double y=0):x(x),y(y){} 31 void in(){inin(shurux),inin(shuruy);x=shurux,y=shuruy;} 32 bool op < (const xl &rhs)const {return x==rhs.x?y<rhs.y:x<rhs.x;} 33 xl op + (const double &rhs)const {return xl(x+rhs,y+rhs);} 34 xl op - (const double &rhs)const {return xl(x-rhs,y-rhs);} 35 xl op * (const double &rhs)const {return xl(x*rhs,y*rhs);} 36 xl op / (const double &rhs)const {return xl(x/rhs,y/rhs);} 37 xl op + (const xl &rhs)const {return xl(x+rhs.x,y+rhs.y);} 38 xl op - (const xl &rhs)const {return xl(x-rhs.x,y-rhs.y);} 39 double angle()const {return atan2(y,x);} 40 double len()const {return sqrt(f(x)+f(y));} 41 }; 42 double D_(const xl &a,const xl &b){return a.x*b.x+a.y*b.y;} 43 double X_(const xl &a,const xl &b){return a.x*b.y-a.y*b.x;} 44 double dis(const xl &a,const xl &b){return sqrt(f(a.x-b.x)+f(a.y-b.y));} 45 double dis2(const xl &a,const xl &b){return (a.x-b.x)*(a.x+b.x)+f(a.y-b.y);} 46 double angle(const xl &a,const xl &b){return acos(D_(a,b)/a.len()/b.len());} 47 struct LINE 48 { 49 xl u,v;double rad; 50 LINE(xl u=xl(0,0),xl v=xl(0,0)):u(u),v(v){} 51 void js(){rad=(v-u).angle();} 52 bool op < (const LINE &rhs)const {return rad<rhs.rad;} 53 }; 54 xl p[100010],ch[100010]; 55 LINE a[100010],b[100010]; 56 bool onleft(const xl &a,const LINE &l){return dcmp(X_(l.v-l.u,a-l.u))>0;} 57 LINE getline2(xl C,xl D) 58 { 59 xl A=p[0],B=p[1]; 60 // if(D_(B-A,D-C)<0)swap(C,D); 61 double a=(-A.y+B.y-C.y+D.y),b=(A.x-B.x+C.x-D.x),c=X_(B,A)+X_(D,C); 62 LINE ret(xl(0,-c/b),xl(1,(-c-a)/b)); 63 if(!dcmp(b))ret=LINE(xl(-c/a,0),xl((-c-b)/a,1)); 64 if(!onleft(p[0],ret)&&!onleft(p[1],ret))swap(ret.u,ret.v); 65 return ret; 66 } 67 LINE getline(xl C,xl D) 68 { 69 LINE l=getline2(C,D); 70 int wocao1=onleft(p[0],l)|onleft(p[1],l); 71 int wocao2=onleft(C,l)|onleft(D,l); 72 if(wocao1^wocao2)return l; 73 else return getline2(D,C); 74 } 75 xl getjd(const LINE &a,const LINE &b) 76 { 77 xl u=a.u-b.u,v1=a.v-a.u,v2=b.v-b.u; 78 double t=X_(v2,u)/X_(v1,v2); 79 return a.u+v1*t; 80 } 81 int n,nn; 82 bool halfplanej(LINE *l) 83 { 84 re(i,0,n-1)l[i].rad=(l[i].v-l[i].u).angle(); 85 sort(l,l+n); 86 int ll=0,rr=0; 87 xl p[n];LINE b[n];b[0]=l[0]; 88 re(i,1,n-1) 89 { 90 while(ll<rr&&!onleft(p[rr-1],l[i]))rr--; 91 while(ll<rr&&!onleft(p[ll],l[i]))ll++; 92 b[++rr]=l[i]; 93 if(!dcmp(X_(b[rr].v-b[rr].u,b[rr-1].v-b[rr-1].u))) 94 { 95 rr--; 96 if(onleft(l[i].u,b[rr]))b[rr]=l[i]; 97 } 98 if(ll<rr)p[rr-1]=getjd(b[rr-1],b[rr]); 99 } 100 while(ll<rr&&!onleft(p[rr-1],b[ll]))rr--; 101 if(rr-ll<=1)return 0; 102 p[rr]=getjd(b[rr],b[ll]); 103 re(i,ll,rr)ch[nn++]=p[i]; 104 ch[nn]=ch[0];return 1; 105 } 106 double zongarea,yaoarea; 107 int main() 108 { 109 inin(n); 110 re(i,0,n-1)p[i].in(); 111 p[n]=p[0]; 112 xl temp=p[1]-p[0]; 113 re(i,1,n-1)zongarea+=abs(X_(temp,p[i+1]-p[i])),temp=temp+(p[i+1]-p[i]); 114 a[0]=LINE(p[0],p[1]); 115 re(i,1,n-1) 116 a[i]=getline2(p[i+1],p[i]); 117 halfplanej(a); 118 temp=ch[1]-ch[0]; 119 re(i,1,nn-1) 120 yaoarea+=abs(X_(temp,ch[i+1]-ch[i])),temp=temp+(ch[i+1]-ch[i]); 121 printf("%.4f",yaoarea/zongarea); 122 return 0; 123 }