【P4196】 [CQOI2006]凸多边形 题解
(参考自《信息学奥赛之数学一本通》)
本题中只要把每条边当作半平面,把每个多边形的每条边做半平面交即可AC。数据范围较小,不必害怕超时。
半平面交的定义和实现方法
直线常用二元一次方程ax+by+c=0表示,我们把平面上的直线及其一侧的部分称为半平面。半平面在平面直角坐标系中可由不等式ax+by+c≥(或≤)0表示。
给定n个形如ax+by+c≤0的半平面,找到所有满足它们的点所组成的点集,称为半平面交。
半平面和半平面交具有以下性质:
(1)半平面或半平面交是一个凸多边形区域,可能无边界,此时增加4个半平面就能保证面积有界;
(2)n个半平面的交是一个至多n条边的凸多边形。
相交后的区域可能是直线、射线、线段或者点,甚至也有可能是空集。
我们先来研究两个凸多边形交的问题。
首先,两个凸多边形的交还是一个凸多边形,一般采用平面扫描法来求两个凸多边形的交。平面扫描法的主要思想是:以两个凸多边形边的交点作为分界点,将边分为内、外两种,内边互相连接,成为所求的多边形。
假设有一条垂直的扫描线,从左到右进行扫描。任何时候,扫描线和两个多边形都至多有4个交点。我们称被扫描到的坐标x为x事件。当然,我们不能扫描所有有理数。但显然,我们能知道下一个x时间将会在当前所扫到的多边形的两条边和垂直扫描线的上下两部分这四条线段的断点和这四条线段的两两之间的交点中选出。整个扫描的时间复杂度为O(n)。
如果求N个半平面的交,则通常使用基于分治思想的D&C算法,分为3个步骤:
第一步,把n个半平面分成两个n/2个半平面的集合;
第二步,把两个子集合递归求解半平面交;
第三步,把前一步算出来的两个半平面交利用平面扫描发求解。
总时间复杂度为O(n log n)。
但是我们还有一种更优的算法:排序增量算法。
先定义一个半平面的极角:比如x-y≥常数的平面,定义它的极角为π/4。
排序增量算法共有以下4个步骤:
(1)把半平面分为两部分,一部分是极角范围之内的(-1/2π,1/2π),另一部分是范围之外的(-π,-1/2π)和(1/2π,π)。在此所有区间都是左开右闭的,下同。
(2)考虑(-1/2π和1/2π)的半平面,把它们按极角排序。对极角相同的半平面,根据常数项保留其中的一个。
(3)从排序后极角最小的两个半平面开始,求出它们的交点并将它们压入一个栈。每次按照极角从小到大的顺序增加一个半平面,算出它与栈顶半平面的交点。如果当前的交点在栈顶两个半平面交点的右边,则让它出栈。注意出栈并不一定只有一次,要继续进行交点检查,如果还有在右边的要继续出栈,直到当前交点在栈顶交点的左边。
(4)相邻半平面的交点组成半个凸多边形,则我们有两个点集:(-1/2π,1/2π)是上半个,(-π,-1/2π)和(1/2π,π)是下半个。
初始的时候,有四个指针p1、p2、p3和p4分别指向上下凸壳的左右端,p1和p3向右走,p2和p4向左走。
在任意时刻,如果最左边的交点不满足p1或p3所在的半平面的限制,这个交点就需要被删除,p1或p3走向它右边的相邻边。我们也可以类似地处理右边的交点,重复运作直到不再有更新出现。
排序增量算法除了第二步的排序之外都是线性的,第二步通常使用快速排序,则总时间复杂度为O(n log n),如果采用基数排序则为O(n)。该算法的常数远小于D&C算法,且编写也更容易,所以一般来说均采用此算法实现。
Code C++
#include<bits/stdc++.h>
using namespace std;
const double eps=1e-13;
struct point
{
double x,y;
point operator -(point &s)
{
return (point){x-s.x,y-s.y};
}
};
double operator *(point a,point b)
{
return a.x*b.y-a.y*b.x;
}
struct line
{
double d;
point a,b;
}l[1005];
bool cmpd(line a,line b)
{
return a.d<b.d;
}
bool bian(point Q,point P1,point P2)
{
return fabs((Q-P1)*(P2-P1)<eps&&min(P1.x,P2.x)-eps<=Q.x&&Q.x-eps<=max(P1.x,P2.x)&&min(P1.y,P2.y)-eps<=Q.y&&Q.y-eps<=max(P1.y,P2.y));
}
point crosp(line a,line b)
{
double s1=(b.a-a.a)*(a.b-a.a),s2=(a.b-a.a)*(b.b-a.a);
return (point){(b.a.x*s2+s1*b.b.x)/(s1+s2),(b.a.y*s2+s1*b.b.y)/(s1+s2)};
}
int n,m,sta[2005];
int main()
{
cin>>n;
for (int i=1;i<=n;i++)
{
int mm;
cin>>mm;
for (int j=1;j<=mm;j++)
{
m++,cin>>l[m].a.x>>l[m].a.y;
l[(j==1)?m+mm-1:m-1].b=l[m].a;
}
}
n=m;
for (int i=1;i<=n;i++)
l[i].d=atan2(l[i].b.y-l[i].a.y,l[i].b.x-l[i].a.x);
sort(&l[1],&l[n+1],cmpd);
for (int i=1;i<=n;i++)
{
for (;sta[0]>=1;sta[0]--)
{
if (fabs(l[i].d-l[sta[sta[0]]].d)<eps)
{
if ((l[sta[sta[0]]].b-l[sta[sta[0]]].a)*(l[i].a-l[sta[sta[0]]].a)<eps) break;
}
else break;
}
for (;sta[0]>=2;sta[0]--)
{
if ((l[i].b-l[i].a)*(crosp(l[sta[sta[0]]],l[sta[sta[0]-1]])-l[i].a)>eps) break;
}
if (fabs(l[i].d-l[sta[sta[0]]].d)>=eps) sta[++sta[0]]=i;
}
int L=1,R=sta[0];
while (L<R)
{
if ((l[sta[L]].b-l[sta[L]].a)*(crosp(l[sta[R]],l[sta[R-1]])-l[sta[L]].a)<eps) R--;
else
{
if ((l[sta[R]].b-l[sta[R]].a)*(crosp(l[sta[L]],l[sta[L+1]])-l[sta[R]].a)<eps) L++;
else break;
}
}
if (R-L<=1)
{
printf("0.000\n");
return 0;
}
double ans=0;
sta[R+1]=sta[L],sta[R+2]=sta[L+1];
for (int i=L;i<=R;i++)
ans+=crosp(l[sta[i]],l[sta[i+1]])*crosp(l[sta[i+1]],l[sta[i+2]])/2;
printf("%.3lf\n",ans);
return 0;
}