[ICPC World Finals 2018][LOJ6409]熊猫保护区(voronoi图)
题面
题解
前置知识
首先转化题意，相当于求多边形内的某一个点，使得这个点到多边形的各定点的距离的最小值最大，输出这个距离。
先求出所有多边形顶点的voronoi图。由voronoi图的性质，对于一个点v和一个多边形顶点u，v到u的距离是v到所有点的距离的最小值当且仅当v在voronoi图中u管辖的区域的内部或边上。所以对于每一个多边形顶点u，“u到voronoi图中它管辖的区域与原多边形的交集中，距离最大的点的距离”的最大值就是答案。
- 性质：对于一个点u和一个多边形A，A中距离u最远的点一定是A的某个顶点。
考虑对于一个多边形顶点u，图形“voronoi图中它管辖的区域与原多边形的交集”的顶点是什么。不可能有除了u以外的其他顶点；可以是voronoi图中的顶点，也可以是voronoi图中的边和多边形边的交点。
所以，原题要求的点只有两种可能性：
- 是voronoi图的某个顶点。
- 是voronoi图中的某条线和原多边形的某条边的交点。
对于1，由于voronoi图只有O(n)个顶点，枚举这些点，再枚举多边形上的点，\(O(n^2)\)暴力即可。
对于2，对于voronoi图中一条线l和多边形某条边的交点p，设l是多边形顶点\(A_i\)与\(A_j\)的中垂线，只用计算p到\(A_i\)或\(A_j\)的距离即可。交点共\(O(n^2)\)个，更新答案用时\(O(1)\)。
综上，本题时间复杂度为\(O(n^2)\)。
代码
#include<bits/stdc++.h>
using namespace std;
#define rg register
#define ld long double
#define In inline
const ld eps = 1e-9;
const ld inf = 1e9;
const int N = 2000;
In int read(){
int s = 0,ww = 1;
char ch = getchar();
while(ch < '0' || ch > '9'){if(ch == '-')ww = -1;ch = getchar();}
while('0' <= ch && ch <= '9'){s = 10 * s + ch - '0';ch = getchar();}
return s * ww;
}
In int sgn(ld x){return x < -eps ? -1 : x > eps;}
In ld sqr(ld x){return x * x;}
struct vec{
ld x,y;
vec(){}
vec(ld _x,ld _y){x = _x,y = _y;}
In friend vec operator + (vec a,vec b){
return vec(a.x + b.x,a.y + b.y);
}
In friend vec operator - (vec a,vec b){
return vec(a.x - b.x,a.y - b.y);
}
In friend vec operator * (vec a,ld k){
return vec(a.x * k,a.y * k);
}
In friend vec operator / (vec a,ld k){
return vec(a.x / k,a.y / k);
}
In friend bool operator == (vec a,vec b){
return sgn(a.x - b.x) == 0 && sgn(a.y - b.y) == 0;
}
In friend ld Dot(vec a,vec b){
return a.x * b.x + a.y * b.y;
}
In friend ld Cross(vec a,vec b){
return a.x * b.y - a.y * b.x;
}
In friend ld len(vec a){
return sqrt(sqr(a.x) + sqr(a.y));
}
In friend bool InUpper(vec a){
int k;
return ((k=sgn(a.y)) > 0 || (k==0&&sgn(a.x)>0));
}
In friend bool prl(vec a,vec b){
return sgn(Cross(a,b)) == 0;
}
In friend bool samedir(vec a,vec b){
return prl(a,b) && sgn(Dot(a,b)) > 0;
}
};
struct line{
vec p,v;
int vx,vy;
line(){}
line(vec _p,vec _v){p = _p,v = _v;}
line(vec _p,vec _v,int _vx,int _vy){p = _p,v = _v,vx = _vx,vy = _vy;}
In friend vec Its(line a,line b){
ld x = Cross(b.v,a.p - b.p),y = Cross(a.v,b.v);
return a.p + a.v * x / y;
}
In friend bool include(line a,vec p){
return sgn(Cross(a.v,p-a.p)) > 0;
}
};
In bool HaveIts(ld l1,ld r1,ld l2,ld r2){
if(l1 > r1)swap(l1,r1);
if(l2 > r2)swap(l2,r2);
if(sgn(r1-l2) < 0 || sgn(r2-l1) < 0)return 0;
return 1;
}
struct seg{
vec p1,p2;
int vx,vy;
seg(){}
seg(vec _p1,vec _p2){p1 = _p1,p2 = _p2;}
seg(vec _p1,vec _p2,int _vx,int _vy){p1 = _p1,p2 = _p2,vx = _vx,vy = _vy;}
In friend bool HaveIts(seg a,seg b){
return HaveIts(a.p1.x,a.p2.x,b.p1.x,b.p2.x)
&& HaveIts(a.p1.y,a.p2.y,b.p1.y,b.p2.y)
&& sgn(Cross(a.p2-a.p1,b.p1-a.p1)) * sgn(Cross(a.p2-a.p1,b.p2-a.p1)) <= 0
&& sgn(Cross(b.p2-b.p1,a.p1-b.p1)) * sgn(Cross(b.p2-b.p1,a.p2-b.p1)) <= 0;
}
In friend vec Its(seg a,seg b){
return Its(line(a.p1,a.p2-a.p1),line(b.p1,b.p2-b.p1));
}
In friend ld dis(vec p,seg a){
return abs(Cross(a.p2-p,a.p1-p) / len(a.p2-a.p1));
}
In friend bool OnSeg(vec p,seg a){
if(p == a.p1 || p == a.p2)return 1;
return sgn(dis(p,a)) == 0 && sgn(Dot(a.p1-p,a.p2-p)) < 0;
}
};
bool inside(vec p,seg l[],int n){//0:out,1:on,2:out
int rt = 0;
for(rg int i = 1;i <= n;i++){
if(OnSeg(p,l[i]))return 1;
seg s = l[i];
if(s.p1.y < s.p2.y)swap(s.p1,s.p2);
if(sgn(p.y-s.p1.y) >= 0 || sgn(p.y-s.p2.y) < 0)continue;
rt ^= (sgn(Cross(p-s.p1,p-s.p2)) > 0);
}
return rt << 1;
}
vec V[N*N+5]; //v图的点，最终是O(n)个，但是过程中可能有很多
seg E[N*N+5];
int Vn,En;
In bool check(line x,line y,line z){
return include(z,Its(x,y));
}
In bool small(line a){
return fabs(a.p.x) < inf / 2 && fabs(a.p.y) < inf / 2;
}
In bool cmp3(line a,line b){
bool k1 = InUpper(a.v),k2 = InUpper(b.v);
if(k1 != k2)return k1 < k2;
return sgn(Cross(a.v,b.v)) > 0;
}
In bool cmp4(line a,line b){
if(samedir(a.v,b.v))return include(b,a.p);
return cmp3(a,b);
}
vec I[N+5];
void HalfPlaneIts(line l[],int ln){
l[++ln] = line(vec(-inf,-inf),vec(1,0));
l[++ln] = line(vec(inf,-inf),vec(0,1));
l[++ln] = line(vec(inf,inf),vec(-1,0));
l[++ln] = line(vec(-inf,inf),vec(0,-1));
sort(l + 1,l + ln + 1,cmp4);
deque<line>q;
q.clear();
for(rg int i = 1;i <= ln;i++){
if(i > 1 && samedir(l[i-1].v,l[i].v))continue;
while(q.size() > 1 && !check(q[q.size()-2],q[q.size()-1],l[i]))q.pop_back();
while(q.size() > 1 && !check(q[1],q[0],l[i]))q.pop_front();
q.push_back(l[i]);
}
while(q.size() > 2 && ! check(q[q.size()-2],q[q.size()-1],q[0]))q.pop_back();
for(rg int i = 0;i < q.size();i++)I[i] = Its(q[i],q[(i+1)%q.size()]);
for(rg int i = 0;i < q.size();i++)if(small(q[i]))E[++En] = seg(I[(i+q.size()-1)%q.size()],I[i],q[i].vx,q[i].vy);
for(rg int i = 0;i < q.size();i++)if(small(q[i]) && small(q[(i+1)%q.size()]))V[++Vn] = I[i];
}
int n;
vec a[N+5]; //多边形的点
seg l[N+5]; //多边形的边
line temp[N+5];
In bool cmp1(vec a,vec b){
int k;
if((k=sgn(a.x-b.x)) != 0)return k < 0;
return sgn(a.y - b.y) < 0;
}
In bool cmp2(vec a,vec b){
return sgn(a.x - b.x) == 0 && sgn(a.y - b.y) == 0;
}
int main(){
n = read();
for(rg int i = 1;i <= n;i++){
int x = read(),y = read();
a[i] = vec(x,y);
}
for(rg int i = 1;i <= n;i++)l[i] = seg(a[i],a[i%n+1]);
for(rg int i = 1;i <= n;i++){
int cnt = 0;
for(rg int j = 1;j <= n;j++)if(i != j)temp[++cnt] = line((a[i]+a[j]) / 2,vec(a[i].y-a[j].y,a[j].x-a[i].x),i,j);
HalfPlaneIts(temp,cnt);
}
sort(V + 1,V + Vn + 1,cmp1);
Vn = unique(V + 1,V + Vn + 1,cmp2) - V - 1;
ld ans = 0;
for(rg int i = 1;i <= Vn;i++){
if(!inside(V[i],l,n))continue;
ld cur = 1e12;
for(rg int j = 1;j <= n;j++)cur = min(cur,len(a[j]-V[i]));
ans = max(ans,cur);
}
for(rg int i = 1;i <= En;i++){
if(InUpper(E[i].p2-E[i].p1))continue;
for(rg int j = 1;j <= n;j++)if(HaveIts(E[i],l[j])){
vec p;
if(prl(E[i].p2-E[i].p1,l[j].p2-l[j].p1))continue;
else p = Its(E[i],l[j]);
ld cur = 1e12;
cur = min(cur,len(a[E[i].vx]-p));
cur = min(cur,len(a[E[i].vy]-p));
ans = max(ans,cur);
}
}
printf("%.9lf\n",(double)ans);
return 0;
}
