【暖*墟】#计算几何# 半平面交的学习与练习

半平面交的相关定义

给出若干个半平面，求它们的公共部分。

其中每个半平面用一条有向直线表示，它的左侧就是它所代表的半平面。

半平面交通常是一个凸多边形，也有时候会得到一个无界多边形，

甚至是线段、直线、点或“空”，无论怎样结果一定是“凸”的。

增量法求半平面交

一般用一个很大的矩形（4个半平面的交）代替“整个平面”，计算结果之后删除这4个人工半平面。

这样，在计算的过程中，每加入一个半平面，都相当于用一条有向直线去切割多边形。

[ 用一条有向直线去切割多边形 ] 的方法：

按逆时针顺序考虑多边形的所有顶点；
保留直线左侧和直线上的点，删除直线右边的点；
如果有向直线和多边形相交时产生了新的点，新点加入新多边形中；
【优化】双端队列：可以动态删除队尾的半平面和队首的半平面。

如图所示，新加入的蓝色线段会删掉凸包首尾，所以要用双端队列。

例题1.【p4196】凸多边形

#include <cmath>
#include <iostream>
#include <cstdio>
#include <string>
#include <cstring>
#include <vector>
#include <algorithm>
#include <stack>
#include <queue>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;

/*【p4196】凸多边形
逆时针给出n个凸多边形的顶点坐标，求它们交的面积。 */

//题意可以转化为：对输入的所有边求半平面交。

//1.计算两个直线的交点（Add函数，用定比分点实现）
//2.计算面积（CalcS函数，用叉积的几何意义实现）
//3.求已有平面与新的一个半平面的交（Cut）

const int N=5019; int n,m,tot,pit_num;

struct Point{ double x,y;
    Point operator - (const Point A){return (Point){x-A.x,y-A.y};}
    Point operator + (const Point A){return (Point){x+A.x,y+A.y};}
    double cross(Point A) {return x*A.y-y*A.x;} //向量叉积
    Point operator * (double t){return (Point){x*t,y*t};} }A[N],B[N],C[N];

void Add(Point a1,Point a2,Point b1,Point b2) //计算两直线交点，放入点坐标数组中
 { Point a=a2-a1,b=b2-b1,c=b1-a1; double t=b.cross(c)/b.cross(a); C[++tot]=a1+a*t; }

void Cut(Point a,Point b){ tot=0; A[pit_num+1]=A[1];
    for(int i=1;i<=pit_num;i++) if((a-A[i]).cross(b-A[i])>=0)
     { C[++tot]=A[i]; if((a-A[i+1]).cross(b-A[i+1])<0) Add(A[i],A[i+1],a,b);}
    else if((a-A[i+1]).cross(b-A[i+1])>=0) Add(A[i],A[i+1],a,b);
    for(int i=1;i<=tot;i++) A[i]=C[i]; pit_num=tot; } //用数组C更新半平面交的答案A

double CalcS() //计算半平面交形成的多边形的总面积
 { double res=0; for(int i=2;i<pit_num;i++) res+=(A[i]-A[1]).cross(A[i+1]-A[1]); return res/2; }

int main(){
    cin>>n>>m; pit_num=m; n--; //n是多边形个数
    for(int i=1;i<=m;i++) cin>>A[i].x>>A[i].y; //初始多边形
    while(n--){ cin>>m>>B[1].x>>B[1].y; B[m+1]=B[1];
        for(int i=2;i<=m;i++) cin>>B[i].x>>B[i].y;
        for(int i=1;i<=m;i++) Cut(B[i],B[i+1]); //计算半平面交
    } printf("%.3f
",CalcS()); return 0; 
}

例题2.【poj2451】Uyuw's Concert

#include <cmath>
#include <iostream>
#include <cstdio>
#include <string>
#include <cstring>
#include <vector>
#include <algorithm>
#include <stack>
#include <queue>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;

/*【poj2451】Uyuw's Concert
给你很多个半平面，求半平面交出来的凸多边形的面积。*/

//题意可以转化为：对输入的所有边求半平面交。

//1.计算两个直线的交点（Add函数，用定比分点实现）
//2.计算面积（CalcS函数，用叉积的几何意义实现）
//3.求已有平面与新的一个半平面的交（Cut）

const int N=500019;

int m,cCnt,curCnt;

struct point{ double x,y; point(double x=0,double y=0):x(x),y(y) {} };

point points[N],p[N],q[N]; //初始化数组points，答案数组p，过程数组q

void getline(point p1,point p2,double &a,double &b,double &c)
 { a=p1.y-p2.y,b=p2.x-p1.x,c=p1.x*p2.y-p2.x*p1.y; }

point intersect(point x,point y,double a,double b,double c){
    double u=fabs(a*x.x+b*x.y+c),v=fabs(a*y.x+b*y.y+c);
    point pt; pt.x=(x.x*v+y.x*u)/(u+v);
    pt.y=(x.y*v+y.y*u)/(u+v); return pt; //求直线的交点并加入点集
}

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]; continue; } //保留
        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; //更新总点数cCnt
}

void solve(){ 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("%.1f
",area); }

int main(){
    points[1]=point(0,0),points[2]=point(10000,0);
    points[3]=point(10000,10000),points[4]=point(0,10000);
    points[5]=p[1]; //初始限定最大矩形
    int n,i;
    while(scanf("%d",&n)!=EOF){
        for(int i=1;i<=4;i++) p[i]=points[i];
        p[5]=p[1]; p[0]=p[4]; cCnt=4;
        for(int i=1;i<=4;i++) //↓↓构造最大矩形并求半平面交
         { double a,b,c; getline(points[i],points[i+1],a,b,c); cut(a,b,c); }
        for(int i=1;i<=n;i++){
            point p1,p2; scanf("%lf%lf%lf%lf",&p1.x,&p1.y,&p2.x,&p2.y);
            double a,b,c; getline(p1,p2,a,b,c); cut(a,b,c);
        } solve(); //求半平面交（凸多边形）的面积
    }
}

——时间划过风的轨迹，那个少年，还在等你