半平面交

    xiaoxiao2021-04-04  188

    转自:http://blog.csdn.net/accry/article/details/6070621

    首先解决问题:什么是半平面? 顾名思义,半平面就是指平面的一半,我们知道,一条直线可以将平面分为两个部分,那么这两个部分就叫做两个半平面。

    然后,半平面怎么表示呢? 二维坐标系下,直线可以表示为ax + by + c = 0,那么两个半平面则可以表示为ax + by + c >= 0 和ax + by + c < 0,这就是半平面的表示方法。 还有,半平面的交是神马玩意? 其实就是一个方程组,让你画出满足若干个式子的坐标系上的区域(类似于线性规划的可行域),方程组就是由类似于上面的这些不等式组成的。

    另外,半平面交可以干什么? 半平面交虽然说是半平面的问题,但它其实就是关于直线的问题。一个一个的半平面其实就是一个一个有方向的直线而已。

    半平面交的一个重要应用就是求多边形的核 。 多边形的核又是神马玩意?  它是平面简单多边形的核是该多边形内部的一个点集,该点集中任意一点与多边形边界上一点的连线都处于这个多边形内部。就是一个在一个房子里面放一个摄像 头,能将所有的地方监视到的放摄像头的地点的集合即为多边形的核。经常会遇到让你判定一个多边形是否有核的问题。

    比如说:

    POJ 3335 Rotating Scoreboard http://acm.pku.edu.cn/JudgeOnline/problem?id=3335

    POJ 1474 Video Surveillance http://acm.pku.edu.cn/JudgeOnline/problem?id=1474

    POJ 1279 Art Gallery http://acm.pku.edu.cn/JudgeOnline/problem?id=1279

    这三个题比较简单,裸的半平面交,下面是核心代码(暨半平面交的模板):

    [cpp]  view plain  copy  print ? /*半平面相交(直线切割多边形)(点标号从1开始)*/   Point points[MAXN],p[MAXN],q[MAXN];   int n;   double r;   int cCnt,curCnt;   inline void getline(Point x,Point y,double &a,double &b,double &c){       a = y.y - x.y;       b = x.x - y.x;       c = y.x * x.y - x.x * y.y;   }   inline void initial(){       for(int i = 1; i <= n; ++i)p[i] = points[i];       p[n+1] = p[1];       p[0] = p[n];       cCnt = n;   }   inline Point intersect(Point x,Point y,double a,double b,double c){       double u = fabs(a * x.x + b * x.y + c);       double v = fabs(a * y.x + b * y.y + c);       return Point( (x.x * v + y.x * u) / (u + v) , (x.y * v + y.y * u) / (u + v) );   }   inline void cut(double a,double b ,double c){       curCnt = 0;       for(int i = 1; i <= cCnt; ++i){           if(a*p[i].x + b*p[i].y + c >= EPS)q[++curCnt] = p[i];           else {               if(a*p[i-1].x + b*p[i-1].y + c > EPS){                   q[++curCnt] = intersect(p[i],p[i-1],a,b,c);               }               if(a*p[i+1].x + b*p[i+1].y + c > EPS){                   q[++curCnt] = intersect(p[i],p[i+1],a,b,c);               }           }       }       for(int i = 1; i <= curCnt; ++i)p[i] = q[i];       p[curCnt+1] = q[1];p[0] = p[curCnt];       cCnt = curCnt;   }   inline void solve(){       //注意:默认点是顺时针,如果题目不是顺时针,规整化方向       initial();       for(int i = 1; i <= n; ++i){           double a,b,c;           getline(points[i],points[i+1],a,b,c);           cut(a,b,c);       }       /*      如果要向内推进r,用该部分代替上个函数      for(int i = 1; i <= n; ++i){          Point ta, tb, tt;          tt.x = points[i+1].y - points[i].y;          tt.y = points[i].x - points[i+1].x;          double k = r / sqrt(tt.x * tt.x + tt.y * tt.y);          tt.x = tt.x * k;          tt.y = tt.y * k;          ta.x = points[i].x + tt.x;          ta.y = points[i].y + tt.y;          tb.x = points[i+1].x + tt.x;          tb.y = points[i+1].y + tt.y;          double a,b,c;          getline(ta,tb,a,b,c);          cut(a,b,c);      }      */       //多边形核的面积       double area = 0;       for(int i = 1; i <= curCnt; ++i)           area += p[i].x * p[i + 1].y - p[i + 1].x * p[i].y;       area = fabs(area / 2.0);       //此时cCnt为最终切割得到的多边形的顶点数,p为存放顶点的数组   }   inline void GuiZhengHua(){        //规整化方向,逆时针变顺时针,顺时针变逆时针       for(int i = 1; i < (n+1)/2; i ++)         swap(points[i], points[n-i]);//头文件加iostream   }   inline void init(){        for(int i = 1; i <= n; ++i)points[i].input();           points[n+1] = points[1];   }  

    半平面交的最终奥义就是求出一个满足条件的凸多边形 ,而解决这个问题的前提就是直线切割多边形 ,让直线不断的去切割当前多边形,然后得到新的多边形,然后继续。。。最终得到一个符合条件的多边形,如果最终的多边形顶点数目为0,那么就说明半平面交无解(半平面交的解可以是一个凸多边形,一条直线,一个点,我们用顶点来记录这些解)。关于直线切割多边形,流程 是这样的: 对 凸多边形(指的是当前多边形)的每一个顶点,如果这个顶点在直线的指定的一侧(暨在该半平面上),那么就将该顶点直接加入到新的多边形中,否则,看与该点 相邻的多边形上的两个点(判断线段是否和直线相交),如果和直线相交,则把交点加入到新的多边形中。这样,就可以得到一个新的凸多边形。关于最初的多边 形,我们可以设一个四方形的区域,比如说(-1000,-1000),(-1000,1000),(1000,1000),(1000,-1000),然 后开始切割。

    POJ 3525 Most Distant Point from the Sea (推荐) http://acm.pku.edu.cn/JudgeOnline/problem?id=3525  求在多边形中画一个圆的最大半径 可以使用半平面交+ 二分法解。对每个距离进行判定(是否存在区域可以放置圆心,多边形的的每条边向内推进

     

    [cpp]  view plain  copy  print ? #include <cstdio>   #include <cstring>   #include <cmath>   #include <iostream>   using namespace std;   #define EPS 1e-10   const int MAXN = 200;   struct Point{       double x,y;       Point(){}       Point(double _x,double _y):x(_x),y(_y){}       void input(){           scanf("%lf%lf",&x,&y);       }   };   Point points[MAXN],p[MAXN],q[MAXN];   int n;   int cCnt,curCnt;   inline void getline(Point x,Point y,double &a,double &b,double &c){        a = y.y - x.y;       b = x.x - y.x;       c = y.x * x.y - x.x * y.y;   }   inline void initial(){       for(int i = 1; i <= n; ++i)p[i] = points[i];       p[n+1] = p[1];       p[0] = p[n];       cCnt = n;   }   inline Point intersect(Point x,Point y,double a,double b,double c){        double u = fabs(a * x.x + b * x.y + c);       double v = fabs(a * y.x + b * y.y + c);       return Point( (x.x * v + y.x * u) / (u + v) , (x.y * v + y.y * u) / (u + v) );   }   inline void cut(double a,double b ,double c){       curCnt = 0;       for(int i = 1; i <= cCnt; ++i){           if(a*p[i].x + b*p[i].y + c >= 0)q[++curCnt] = p[i];           else {               if(a*p[i-1].x + b*p[i-1].y + c > 0){                   q[++curCnt] = intersect(p[i],p[i-1],a,b,c);               }               if(a*p[i+1].x + b*p[i+1].y + c > 0){                   q[++curCnt] = intersect(p[i],p[i+1],a,b,c);               }           }       }       for(int i = 1; i <= curCnt; ++i)p[i] = q[i];       p[curCnt+1] = q[1];p[0] = p[curCnt];       cCnt = curCnt;   }   inline int solve(double r){       initial();       for(int i = 1; i <= n; ++i){           Point ta, tb, tt;           tt.x = points[i+1].y - points[i].y;           tt.y = points[i].x - points[i+1].x;           double k = r / sqrt(tt.x * tt.x + tt.y * tt.y);           tt.x = tt.x * k;           tt.y = tt.y * k;           ta.x = points[i].x + tt.x;           ta.y = points[i].y + tt.y;           tb.x = points[i+1].x + tt.x;           tb.y = points[i+1].y + tt.y;           double a,b,c;           getline(ta,tb,a,b,c);           cut(a,b,c);       }       if(cCnt <= 0)return 0;       return 1;   }   int main(){       while(scanf("%d",&n) != EOF){           if(n == 0)break;           for(int i = 1; i <= n; ++i)points[i].input();          for(int i = 1; i < (n+1)/2; i ++)               swap(points[i], points[n-i]);           points[n+1] = points[1];           double high = 10000000;           double low = 0,mid;           while(low + EPS <= high){               mid = (high + low)/2.0;               if(solve(mid))low = mid;               else high = mid;           }           printf("%lf/n",high);       }       return 0;   }  

     

    POJ 3384 Feng Shui (推荐) http://acm.pku.edu.cn/JudgeOnline/problem?id=3384  半平面交实际应用,用两个圆覆盖一个多边形,问最多能覆盖多边形的面积。 解法:用半平面交将多边形的每条边一起向“内”推进R,得到新的多边形,然后求多边形的最远两点

     

    [cpp]  view plain  copy  print ? #include <cstdio>   #include <cstring>   #include <cmath>   #include <iostream>   using namespace std;   #define EPS 1e-10   const int MAXN = 300;   struct Point{       double x,y;       Point(){}       Point(double _x,double _y):x(_x),y(_y){}       void input(){           scanf("%lf%lf",&x,&y);       }   };   Point points[MAXN],p[MAXN],q[MAXN];   int n;   double r;   int cCnt,curCnt;   void getline(Point x,Point y,double &a,double &b,double &c){    a = y.y - x.y;    b = x.x - y.x;    c = y.x * x.y - x.x * y.y;   }   inline void initial(){       for(int i = 1; i <= n; ++i)p[i] = points[i];       p[n+1] = p[1];       p[0] = p[n];       cCnt = n;   }   inline Point intersect(Point x,Point y,double a,double b,double c){        double u = fabs(a * x.x + b * x.y + c);       double v = fabs(a * y.x + b * y.y + c);       return Point( (x.x * v + y.x * u) / (u + v) , (x.y * v + y.y * u) / (u + v) );   }   inline void cut(double a,double b ,double c){       curCnt = 0;       for(int i = 1; i <= cCnt; ++i){           if(a*p[i].x + b*p[i].y + c >= 0)q[++curCnt] = p[i];           else {               if(a*p[i-1].x + b*p[i-1].y + c > 0){                   q[++curCnt] = intersect(p[i],p[i-1],a,b,c);               }               if(a*p[i+1].x + b*p[i+1].y + c > 0){                   q[++curCnt] = intersect(p[i],p[i+1],a,b,c);               }           }       }       for(int i = 1; i <= curCnt; ++i)p[i] = q[i];       p[curCnt+1] = q[1];p[0] = p[curCnt];       cCnt = curCnt;   }   inline double pdis(Point a,Point b){       return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));   }   int main(){       scanf("%d%lf",&n,&r);       for(int i = 1; i <= n; ++i)points[i].input();       points[n+1] = points[1];       initial();       for(int i = 1; i <= n; ++i){           Point ta, tb, tt;           tt.x = points[i+1].y - points[i].y;           tt.y = points[i].x - points[i+1].x;           double k = r / sqrt(tt.x * tt.x + tt.y * tt.y);           tt.x = tt.x * k;           tt.y = tt.y * k;           ta.x = points[i].x + tt.x;           ta.y = points[i].y + tt.y;           tb.x = points[i+1].x + tt.x;           tb.y = points[i+1].y + tt.y;           double a,b,c;           getline(ta,tb,a,b,c);           cut(a,b,c);       }       int ansx = 0,ansy = 0;       double res = 0;       for(int i = 1; i <= cCnt; ++i){           for(int j = i + 1; j <= cCnt; ++j){               double tmp = pdis(p[i],p[j]);               if(tmp > res){                   res = tmp;                   ansx = i;                   ansy = j;               }           }       }       printf("%.4lf %.4lf %.4lf %.4lf/n",p[ansx].x,p[ansx].y,p[ansy].x,p[ansy].y);       return 0;   }  

     

    POJ 1755 Triathlon (推荐) http://acm.pku.edu.cn/JudgeOnline/problem?id=1755  半平面交判断不等式是否有解。注意特殊情况

     

    [cpp]  view plain  copy  print ? #include <cstdio>   #include <cstring>   #include <cmath>   using namespace std;   #define EPS 1e-8   const int MAXN = 110;   struct Point{       double x,y;       Point(){}       Point(double _x,double _y):x(_x),y(_y){}   };   struct Node{       double x,y,z;       void input(){           scanf("%lf%lf%lf",&x,&y,&z);       }   };   inline void getabc(double &a,double &b,double &c,Node x1,Node x2){       a = 1.0/x2.x - 1.0/x1.x ;       b = 1.0/x2.y - 1.0/x1.y;       c = 1.0/x2.z - 1.0/x1.z;   }   Point points[MAXN],p[MAXN],q[MAXN];   int n;   int cCnt,curCnt;   int flag ;   int sig(double k){    return (k < -EPS) ? -1 : (k > EPS);   }   inline void getline(Point x,Point y,double &a,double &b,double &c){       a = y.y - x.y;       b = x.x - y.x;       c = y.x * x.y - x.x * y.y;   }   inline void initial(){       p[1] = Point(0,0);       p[2] = Point(0,1000000);       p[3] = Point(1000000,1000000);       p[4] = Point(1000000,0);       p[5] = p[1];       p[0] = p[4];       cCnt = 4;   }   inline Point intersect(Point x,Point y,double a,double b,double c){       double u = fabs(a * x.x + b * x.y + c);       double v = fabs(a * y.x + b * y.y + c);       return Point( (x.x * v + y.x * u) / (u + v) , (x.y * v + y.y * u) / (u + v) );   }   inline void cut(double a,double b ,double c){       curCnt = 0;       if(a == 0 && b == 0 && c == 0){           flag = 0;           return ;       }       for(int i = 1; i <= cCnt; ++i){           if(sig(a*p[i].x + b*p[i].y + c) >= 0)q[++curCnt] = p[i];           //注意这里一定是大于0,而不是大于等于0,直接排除掉a = b = c = 0这种情况           else {               if(sig(a*p[i-1].x + b*p[i-1].y + c) > 0){                   q[++curCnt] = intersect(p[i],p[i-1],a,b,c);               }               if(sig(a*p[i+1].x + b*p[i+1].y + c) > 0){                   q[++curCnt] = intersect(p[i],p[i+1],a,b,c);               }           }       }       double area = 0;       for(int i = 1; i <= curCnt; ++i)           area += p[i].x * p[i + 1].y - p[i + 1].x * p[i].y;       area = fabs(area / 2.0);       if(area <= EPS)flag = 0;       for(int i = 1; i <= curCnt; ++i)p[i] = q[i];       p[curCnt+1] = q[1];p[0] = p[curCnt];       cCnt = curCnt;   }   Node nodes[MAXN];   int N;   int main(){       scanf("%d",&N);       for(int i = 1; i <= N; ++i)           nodes[i].input();       bool ok;       double a,b,c;       for(int i = 1; i <= N; ++i){           ok = 1;           initial();           flag = 1;           for(int j = 1; j <= N; j++){               if(i == j)continue;               getabc(a,b,c,nodes[i],nodes[j]);               cut(a,b,c);               if(flag == 0){                   ok = 0;                   break;               }           }           if(ok)printf("Yes/n");           else printf("No/n");       }       return 0;   }  

     

     

    POJ 2540 Hotter Colder http://acm.pku.edu.cn/JudgeOnline/problem?id=2540 半平面交求线性规划可行区域的面积,需要求两点的中垂线。 注意直线切割多边形的时候,向新的点数组(q)里加点的时候,点的顺序必须和原来保持一致(顺时针)

     

     

    [cpp]  view plain  copy  print ? #include <cstdio>   #include <cstring>   #include <cmath>   using namespace std;   #define EPS 1e-6   const int MAXN = 200;   struct Point{       double x,y;       Point(){}       Point(double _x,double _y):x(_x),y(_y){}       void input(){           scanf("%lf%lf",&x,&y);       }   };   Point p[MAXN],q[MAXN];   int cCnt ,curCnt;   inline Point intersect(Point x,Point y,double a,double b,double c){        double u = fabs(a * x.x + b * x.y + c);       double v = fabs(a * y.x + b * y.y + c);       return Point( (x.x * v + y.x * u) / (u + v) , (x.y * v + y.y * u) / (u + v) );   }   int sig(double k){    return (k < -EPS) ? -1 : (k > EPS);   }   inline void gethalf(Point p1,Point p2,double &a,double &b,double &c,                        int flag){       int fhalf;       a = p2.x - p1.x;       b = p2.y - p1.y;       c = (p1.y*p1.y - p2.y*p2.y - p2.x*p2.x + p1.x*p1.x)/2;       double tmp = a*p2.x + b*p2.y + c;       if(flag == 1){           if(sig(tmp) >= 0)fhalf = 1;           else fhalf = -1;       }       if(flag == -1){           if(sig(tmp) >= 0)fhalf = -1;           else fhalf = 1;       }       if(fhalf == -1){           a = -a;           b = -b;           c = -c;       }   }   inline void cut(double a,double b,double c){       curCnt = 0;           for(int i = 1; i <= cCnt; ++i){               if(sig(a*p[i].x + b*p[i].y + c) >= 0)q[++curCnt] = p[i];               else {                   if(sig(a*p[i-1].x + b*p[i-1].y + c) > 0){                       q[++curCnt] = intersect(p[i-1],p[i],a,b,c);                   }                   if(sig(a*p[i+1].x + b*p[i+1].y + c) > 0){                       q[++curCnt] = intersect(p[i+1],p[i],a,b,c);                   }               }           }       for(int i = 1; i <= curCnt; ++i)           p[i] = q[i];       p[curCnt + 1] = p[1];       p[0] = p[curCnt];       cCnt = curCnt;   }   int main(){       p[1] = Point(0,0);       p[2] = Point(10,0);       p[3] = Point(10,10);       p[4] = Point(0,10);       p[5] = p[1];       p[0] = p[4];       cCnt = 4;       double a1,b1;       char op[20];       bool mark = 1;       Point last = Point(0,0);       while(scanf("%lf%lf%s",&a1,&b1,&op) != EOF){           Point current = Point(a1,b1);           if(!mark){               puts("0.00");               continue;           }           if(!strcmp(op,"Same")){               puts("0.00");               mark = 0;               continue;           }           int fhalf = 0,flag = 0;           if(!strcmp(op,"Hotter"))flag = 1;           else if(!strcmp(op,"Colder"))flag = -1;           double a,b,c;           gethalf(last,current,a,b,c,flag);           cut(a,b,c);           double area = 0;           for(int i = 1; i <= cCnt; ++i)               area += p[i].x * p[i + 1].y - p[i + 1].x * p[i].y;           area = fabs(area / 2.0);           printf("%.2lf/n",area);           last = current;       }       return 0;   }  

     

     

     

    另外还有一个数据比较多的,必须用nlogn的算法解决的半平面交的问题:

    http://acm.pku.edu.cn/JudgeOnline/problem?id=2451

    Zzy 专为他那篇nlogn算法解决半平面交问题的论文而出的题目,至今不会自己写。

    关于求多边形内核的算法

    什么是多边形的内核?

    它是平面简单多边形的核是该多边形内部的一个点集,该点集中任意一点与多边形边界上一点的连线都处于这个多边形内部。就是一个在一个房子里面放一个摄像 头,能将所有的地方监视到的放摄像头的地点的集合即为多边形的核。

     

             

    如上图,第一个图是有内核的,比如那个黑点,而第二个图就不存在内核了,无论点在哪里,总有地区是看不到的。

     

    那么,如何求得这个内核区间呢?通常的算法是用两点的直线去不断切割多边形,切割到最后剩下的,就是内核区间了。

    我们都知道一条直线可以将平面切割成两个区域,假设直线方程为

    ax+by+c==0,那么,两个平面可分别表示成ax+by+c>=0 和 ax+by+c<0

     

    具体如何用程序实现直线对多边形的切割呢?

    流程是这样的:

    1、          用一个顺时针或者逆时针的顺序,将最初的多边形的点集储存起来。

    2、          按顺序取连续的两个点组成一条直线,用这条直线来切割原先的多边形

    我首先假设点是顺时针储存的,如图:

    此时,多边形的点集是{1,2,3,4,5,6,7,8,9,10}

     

    取点1,和点2组成直线ax+by+c==0,这时候,将点集中的点一次带入方程ax+by+c,得到的值都将会是大于等于0的,说明所有的点都在该直线的同一侧,继续保持点集不变

     

    取点2和点3组成直线,同样,将点集中的点依次带入方程ax+by+c中,此时,4和5两个点的结果是小于0的,而其他的点的值依旧是大于等于0,这时候说明4和5两个点被切割出了该多边形,于是现在点集只剩下{1,2,3,6,7,8,9,10,X},(X是直线23和直线56的交点)

     

    依次类推,一直执行到点10和点1,那么内核的集合就得到了。

     

    值得说明的是,这个例子的图形比较特殊,全是直角,如果图形比较随意,那么当某一个点被断定在多边形区间之外的时候,我们还应该考虑它和它相邻的两个点各自组成的直线和ax+by+c有没有交点,有交点的话,更新的点集中还应该加上这些交点,比如例子中执行完点2和点3组成的直线后,点集是{1,2,3,6,7,8,X},其中3和X就是这样的结果

     

    还有,为什么将所有的点依次执行一遍,然后取剩下的某一边的点构成新的点集就够了呢?答案是,点是顺时针或者逆时针给出的~~~

    转载请注明原文地址: https://ju.6miu.com/read-666159.html

    最新回复(0)