ECNU 1624 求交集多边形面积

其实我是想做hdu3060的。。。写完了发现。。。这题坑姐，还有凹的情况。。。无奈先找个OJ把俩凸的代码给交了。。

求俩凸多边形面积。水题啊。本来我想了一种方法，想着能水过去呢。一直WA。后来用半平面交算了。。。好久没写半平面交的题了。。。

刚改了点BUG，我的那个方法也水过去啦！！！(*^__^*) 嘻嘻……

我的是M*N的算法，半平面交是log(N+M)*(N+M)的。

我这个就很好想啦，既然是求面积的交，因为是凸多边形，所以两个的交集一定也是凸的（如果存在的话）。

而且稍微想一下就很清楚，交得的凸多边形的顶点要么是交点，要么是在另一个凸多边形里的点。所以把这些个点弄到一个集合里，极角排序，然后求面积即可~~~
需要特判是否一个完全包含另一个，还有不想交的情况~

中间一直WA是因为，求得的点集如果有重点，那么用极角排序（用atan2）是很不靠谱的！！！因为两个相同点的atan2的结果是0！而极角排序是从负值到正值，如果俩点相同，经过极角排序后，它俩不在一起了，哭。所以，在极角排序前去重下就没问题啦。

有很大改进空间滴。。

我的N*M做法

#include <set>#include <map>#include <queue>#include <stack>#include <math.h>#include <stdio.h>#include <stdlib.h>#include <iostream>#include <limits.h>#include <string.h>#include <string>#include <algorithm>#define MID(x,y) ( ( x + y ) >> 1 )#define L(x) ( x << 1 )#define R(x) ( x << 1 | 1 )#define FOR(i,s,t) for(int i=(s); i<(t); i++)#define BUG puts("here!!!")#define STOP system("pause")#define file_r(x) freopen(x, "r", stdin)#define file_w(x) freopen(x, "w", stdout)using namespace std;const int MAX = 510;const double eps = 1e-6;bool dy(double x,double y){return x > y + eps;}// x > y bool xy(double x,double y){return x < y - eps;}// x < y bool dyd(double x,double y){ return x > y - eps;}// x >= y bool xyd(double x,double y){return x < y + eps;} // x <= y bool dd(double x,double y) {return fabs( x - y ) < eps;}  // x == ystruct point{double x, y;void get(){scanf("%lf%lf", &x, &y);}bool operator==(const point &a){return dd(a.x, x) && dd(a.y, y);}};double crossProduct(point a,point b,point c)//向量 ac 在 ab 的方向 顺时针是正 {return (c.x - a.x)*(b.y - a.y) - (b.x - a.x)*(c.y - a.y);}point l2l_inst_p(point u1,point u2,point v1,point v2){point ans = u1;double t = ((u1.x - v1.x)*(v1.y - v2.y) - (u1.y - v1.y)*(v1.x - v2.x))/((u1.x - u2.x)*(v1.y - v2.y) - (u1.y - u2.y)*(v1.x - v2.x));ans.x += (u2.x - u1.x)*t;ans.y += (u2.y - u1.y)*t;return ans;}bool onSegment(point a, point b, point c){if( dd(crossProduct(a,b,c),0.0) && dyd(c.x,min(a.x,b.x)) && xyd(c.x,max(a.x,b.x)) && dyd(c.y,min(a.y,b.y)) && xyd(c.y,max(a.y,b.y)) )return true;return false;}bool s2s_inst(point p1,point p2, point p3, point p4) {double d1 = crossProduct(p3,p4,p1);double d2 = crossProduct(p3,p4,p2);double d3 = crossProduct(p1,p2,p3);double d4 = crossProduct(p1,p2,p4);if( xy(d1 * d2, 0.0) && xy(d3 * d4, 0.0) ) return true;return false;}point p1[MAX], p2[MAX];bool pin_convexh(point *p,int n,point a){p[n] = p[0]; p[n+1] = p[1];for(int i=0; i<n; i++)if( xy(crossProduct(p[i],p[i+1],a)*crossProduct(p[i+1],p[i+2],a),0.0) )return false;return true;} double area_polygon(point p[],int n){if( n < 3 ) return 0.0; double s = 0.0;for(int i=0; i<n; i++)s += p[(i+1)%n].y * p[i].x - p[(i+1)%n].x * p[i].y;return fabs(s)/2.0;}bool ainb(point *a, int n, point *b, int m){FOR(i, 0, n)if( !pin_convexh(b, m, a[i]) )return false;return true;}bool inst(point *a, int n, point *b, int m){FOR(i, 0, n)FOR(k, 0, m)if( s2s_inst(a[i], a[i+1], b[k], b[k+1]) )return true;return false;}point t[MAX*2], C;bool cmp(point a,point b)  {  double t1 = atan2(a.y - C.y, a.x - C.x);double t2 = atan2(b.y - C.y, b.x - C.x);    if( dd(t1, t2) ) return xy(fabs(a.x),fabs(b.x));    return xy(t1, t2);  }bool cmp_p(point a, point b){if( dd(a.x, b.x) )return xy(a.y, b.y);return xy(a.x, b.x);}bool cmp_equal(point a, point b){return a == b;}double solve(int n, int m){double area1 = area_polygon(p1, n);double area2 = area_polygon(p2, m);if( ainb(p1, n, p2, m) || ainb(p2, m, p1, n) )// 如果一个在另一个里面，答案是面积最小的 return min(area1, area2);if( !inst(p1, n, p2, m) )// 如果两个不相交return 0;// 如果两个相交，纠结啊int cnt = 0;FOR(i, 0, n)FOR(k, 0, m)if( s2s_inst(p1[i], p1[i+1], p2[k], p2[k+1]) )t[cnt++] = l2l_inst_p(p1[i], p1[i+1], p2[k], p2[k+1]);FOR(i, 0, n)if( pin_convexh(p2, m, p1[i]) )t[cnt++] = p1[i];FOR(i, 0, m)if( pin_convexh(p1, n, p2[i]) )t[cnt++] = p2[i];sort(t, t+cnt, cmp_p);cnt = unique(t, t+cnt, cmp_equal) - t;C = t[0];sort(t+1, t+cnt, cmp);cnt = unique(t, t+cnt, cmp_equal) - t;double area = area_polygon(t, cnt);return area;} int main(){int n, m;while( ~scanf("%d", &n) ){FOR(i, 0, n)p1[i].get();scanf("%d", &m);FOR(i, 0, m)p2[i].get();p1[n] = p1[0]; p2[m] = p2[0];double ans = solve(n, m);printf("%.2lf\n", ans);}return 0;}

半平面交做法

#include <set>#include <map>#include <queue>#include <stack>#include <math.h>#include <stdio.h>#include <stdlib.h>#include <iostream>#include <limits.h>#include <string.h>#include <string>#include <algorithm>#define MID(x,y) ( ( x + y ) >> 1 )#define L(x) ( x << 1 )#define R(x) ( x << 1 | 1 )#define FOR(i,s,t) for(int i=(s); i<(t); i++)#define BUG puts("here!!!")#define STOP system("pause")#define file_r(x) freopen(x, "r", stdin)#define file_w(x) freopen(x, "w", stdout)using namespace std;const int MAX = 510;const double eps = 1e-8;bool dy(double x,double y){return x > y + eps;}// x > y bool xy(double x,double y){return x < y - eps;}// x < y bool dyd(double x,double y){ return x > y - eps;}// x >= y bool xyd(double x,double y){return x < y + eps;} // x <= y bool dd(double x,double y) {return fabs( x - y ) < eps;}  // x == ystruct point{double x, y;void get(){scanf("%lf%lf", &x, &y);}bool operator==(const point &a){return dd(a.x, x) && dd(a.y, y);}};double crossProduct(point a,point b,point c)//向量 ac 在 ab 的方向 顺时针是正 {return (c.x - a.x)*(b.y - a.y) - (b.x - a.x)*(c.y - a.y);}double area_polygon(point p[],int n){if( n < 3 ) return 0.0; double s = 0.0;for(int i=0; i<n; i++)s += p[(i+1)%n].y * p[i].x - p[(i+1)%n].x * p[i].y;return fabs(s)/2.0;}point t[MAX*2], p1[MAX], p2[MAX];struct line{point a, b;double ang;};line ln[MAX*2], deq[MAX*2];bool cmp_equal(point a, point b){return a == b;}bool equal_ang(line a,line b)// 第一次unique的比较函数 {return dd(a.ang, b.ang);}bool cmphp(line a,line b)// 排序的比较函数 {if( dd(a.ang,b.ang) ) return xy(crossProduct(b.a,b.b,a.a),0.0);return xy(a.ang,b.ang);}bool equal_p(point a,point b)//第二次unique的比较函数 {return dd(a.x,b.x) && dd(a.y,b.y);}void makeline_hp(double x1,double y1,double x2,double y2,line &l){l.a.x = x1; l.a.y = y1; l.b.x = x2; l.b.y = y2;l.ang = atan2(y2 - y1,x2 - x1);}void makeline_hp(point a,point b,line &l) // 线段(向量ab)左侧侧区域有效 {l.a = a; l.b = b;l.ang = atan2(b.y - a.y,b.x - a.x);// 如果是右侧区域，改成a.y - b.y,a.x - b.x }bool parallel(line u,line v){return dd( (u.a.x - u.b.x)*(v.a.y - v.b.y) - (v.a.x - v.b.x)*(u.a.y - u.b.y) , 0.0 );}point l2l_inst_p(line l1,line l2){point ans = l1.a;double t = ((l1.a.x - l2.a.x)*(l2.a.y - l2.b.y) - (l1.a.y - l2.a.y)*(l2.a.x - l2.b.x))/   ((l1.a.x - l1.b.x)*(l2.a.y - l2.b.y) - (l1.a.y - l1.b.y)*(l2.a.x - l2.b.x));ans.x += (l1.b.x - l1.a.x)*t;ans.y += (l1.b.y - l1.a.y)*t;return ans;}void inst_hp_nlogn(line *ln,int n,point *s,int &len){len = 0;sort(ln,ln+n,cmphp);n = unique(ln,ln+n,equal_ang) - ln;int bot = 0,top = 1;deq[0] = ln[0]; deq[1] = ln[1];for(int i=2; i<n; i++){if( parallel(deq[top],deq[top-1]) || parallel(deq[bot],deq[bot+1]) )return ;while( bot < top && dy(crossProduct(ln[i].a,ln[i].b,l2l_inst_p(deq[top],deq[top-1])),0.0) )top--;while( bot < top && dy(crossProduct(ln[i].a,ln[i].b,l2l_inst_p(deq[bot],deq[bot+1])),0.0) )bot++;deq[++top] = ln[i];}while( bot < top && dy(crossProduct(deq[bot].a,deq[bot].b,l2l_inst_p(deq[top],deq[top-1])),0.0) )top--;while( bot < top && dy(crossProduct(deq[top].a,deq[top].b,l2l_inst_p(deq[bot],deq[bot+1])),0.0) )bot++;if( top <= bot + 1 ) return ;for(int i=bot; i<top; i++)s[len++] = l2l_inst_p(deq[i],deq[i+1]);if( bot < top + 1 ) s[len++] = l2l_inst_p(deq[bot],deq[top]);len = unique(s,s+len,equal_p) - s;} double solve(int n, int m){int cnt = 0;FOR(i, 0, n)makeline_hp(p1[i], p1[i+1], ln[cnt++]);FOR(i, 0, m)makeline_hp(p2[i], p2[i+1], ln[cnt++]);int len;inst_hp_nlogn(ln, cnt, t, len);if( len == 0 ) return 0;double area = area_polygon(t, len);return area;} int main(){int n, m;while( ~scanf("%d", &n) ){for(int i=n-1; i>=0; i--)p1[i].get();scanf("%d", &m);for(int i=m-1; i>=0; i--)p2[i].get();p1[n] = p1[0]; p2[m] = p2[0];double ans = solve(n, m);printf("%.2lf\n", ans);}return 0;}