# 引入Voronoi图，定义法

//n=2000,time=1800ms
#include<cstdio>
#include<cmath>
#include<iostream>
#include<iomanip>
#include<algorithm>
#define double long double
#define TEST
using namespace std;
const int N=2009;
const double PI=acos(-1),EPS=1e-10;
struct Point{
double x,y;
Point(){x=y=0;}
Point(double a){x=a;y=0;}
Point(double a,double b){x=a;y=b;}
friend istream&operator>>(istream&is,Point&a){return is>>a.x>>a.y;}
friend ostream&operator<<(ostream&os,Point a){return os<<a.x<<','<<a.y;}
Point operator-(){return Point(-x,-y);}
Point operator+(Point a){return Point(x+a.x,y+a.y);}
Point operator-(Point a){return Point(x-a.x,y-a.y);}
Point operator*(double a){return Point(x*a,y*a);}
Point operator/(double a){return Point(x/a,y/a);}
friend Point operator*(double a,Point b){b.x*=a,b.y*=a;return b;}
Point&operator+=(Point a){x+=a.x,y+=a.y;return*this;}
Point&operator-=(Point a){x-=a.x,y-=a.y;return*this;}
Point&operator*=(double a){x*=a,y*=a;return*this;}
Point&operator/=(double a){x/=a,y/=a;return*this;}
double operator*(Point a){return x*a.x+y*a.y;}
double operator%(Point a){return x*a.y-y*a.x;}
double Len(){return sqrt(x*x+y*y);}
friend double Dis(Point a,Point b){return (a-b).Len();}
Point Turn90(){return Point(-y,x);}
}a[N];
struct Line{
Point p,v;double ang;
Line(){}
Line(Point a,Point b){p=a,v=b-a,ang=atan2(v.y,v.x);}
friend ostream&operator<<(ostream&os,Line a){return os<<a.p<<' '<<a.v;}
bool operator<(Line&a){return ang<a.ang;}
bool Right(Point&a){return v%(a-p)<EPS;}
friend Point Cross(Line a,Line b){return a.p+a.v*(b.v%(b.p-a.p))/(b.v%a.v);}
}l[N];
int n,w,h;
Line q[N],*qq;
Point k[N],*kk;
int HalfPlane(Line*a,int n){//半平面交
int h=0,t=0;
sort(a,a+n);
q[0]=a[0];
for(int i=1;i<n;++i){
while(h<t&&a[i].Right(k[t-1]))--t;
while(h<t&&a[i].Right(k[h]))++h;
if(fabs(q[t].ang-a[i].ang)>EPS)q[++t]=a[i];
else if(a[i].Right(q[t].p))q[t]=a[i];
if(h<t)k[t-1]=Cross(q[t-1],q[t]);
}
while(h<t&&q[h].Right(k[t-1]))--t;
k[t]=Cross(q[h],q[t]);
qq=q+h;
kk=k+h;
return t-h+1;
}
int main(){
freopen("in.in","r",stdin);
freopen("halfplane.out","w",stdout);
cin>>n>>w>>h;
for(int i=1;i<=n;++i)
cin>>a[i];
double ans=0;
for(int i=1;i<=n;++i){
int m=0;
for(int j=1;j<=n;++j)
if(i!=j){
Point mid=(a[i]+a[j])/2,v=a[j]-a[i];
l[m++]=Line(mid,mid+v.Turn90());
}
l[m++]=Line(Point(0,0),Point(w,0));
l[m++]=Line(Point(w,0),Point(w,h));
l[m++]=Line(Point(w,h),Point(0,h));
l[m++]=Line(Point(0,h),Point(0,0));
m=HalfPlane(l,m);
#ifdef TEST//测试生成多边形的边数，检查是否成功处理退化情况
cout<<m<<endl;
#endif
for(int j=0;j<m;++j)
ans=max(ans,Dis(a[i],kk[j]));
}
cout<<ans<<endl;
return 0;
}


# 引入Delaunay三角网，逐点插入法

Voronoi图是平面图，Delaunay三角网与它互为对偶图（对于Voronoi图的每条边，连接其相邻两个泰森多边形的基点）

（暂不讨论特殊情况：多个基点共圆，Voronoi图的每个顶点与更多边相连，多个Delaunay三角形外接圆圆心重合）

主函数
初始化一个极大三角形
枚举点i：1~n
枚举三角形，找到i所在的三角形abc
加入三角形abi，bci，aci
FlipTest(a,b,i)
FlipTest(b,c,i)
FlipTest(a,c,i)

找到ab边所属另一侧的三角形的第三个顶点d
如果i在abd的外接圆内（不满足空圆性）
删除abi，abd
加入adi，bdi（形象理解，相当于把四边形adbi中间横着的ab边flip了一下）
FlipTest(a,d,i)
FlipTest(b,d,i)


## 不依赖Delaunay三角网的逐点插入法

初始化极大边界

枚举点j：1~i-1
找到与点i最近的点near
枚举点j：near开始，再次到near退出
作ij中垂线，与面j交于点p1,p2
更新面i（加入边p2p1）
更新面j（加入边p1p2，删除外侧多余的边）
j:=p2所在边的另一侧面对应的点


update: 基本解决了，注意的地方见注释

//n=2000,time=120ms
//n=10000,time=2000ms
#include<cstdio>
#include<cmath>
#include<iostream>
#include<algorithm>
#define double long double
#define TEST
using namespace std;
const int N=10009,E=30*N;
const double PI=acos(-1),EPS=1e-9,INF=1e4*(1+1/sqrt(2));
struct Point{
double x,y;
Point(){x=y=0;}
Point(double a){x=a;y=0;}
Point(double a,double b){x=a;y=b;}
friend istream&operator>>(istream&is,Point&a){return is>>a.x>>a.y;}
friend ostream&operator<<(ostream&os,Point a){return os<<a.x<<','<<a.y;}
Point operator-(){return Point(-x,-y);}
Point operator+(Point a){return Point(x+a.x,y+a.y);}
Point operator-(Point a){return Point(x-a.x,y-a.y);}
Point operator*(double a){return Point(x*a,y*a);}
Point operator/(double a){return Point(x/a,y/a);}
friend Point operator*(double a,Point b){b.x*=a,b.y*=a;return b;}
Point&operator+=(Point a){x+=a.x,y+=a.y;return*this;}
Point&operator-=(Point a){x-=a.x,y-=a.y;return*this;}
Point&operator*=(double a){x*=a,y*=a;return*this;}
Point&operator/=(double a){x/=a,y/=a;return*this;}
double operator*(Point a){return x*a.x+y*a.y;}
double operator%(Point a){return x*a.y-y*a.x;}
double Len(){return sqrt(x*x+y*y);}
friend double Dis(Point a,Point b){return (a-b).Len();}
Point Turn(double r){double c=cos(r),s=sin(r);return Point(x*c-y*s,y*c+x*s);}
Point Turn90(){return Point(-y,x);}
bool operator==(Point a){return fabs(x-a.x)<EPS&&fabs(y-a.y)<EPS;}
bool operator!=(Point a){return fabs(x-a.x)>EPS||fabs(y-a.y)>EPS;}
};
int ecnt=1,ehead[N],face[E],suc[E];
Point a[N],p[E];

Point Intersect(Point a,Point av,Point b,Point bv){//直线求交点
return a+av*(bv%(b-a))/(bv%av);
}
void AddEdge(Point a,Point b,int af,int bf){//连边
p[++ecnt]=a;face[ecnt]=af;
p[++ecnt]=b;face[ecnt]=bf;
}
int Cut(Point m,Point mv,int e,int f){//有向直线(m,mv)切割边e所在面，另一侧划分给面f，返回翻边
int ein=e,eout,esuc;//切入边，切出边，后继
Point pin,pout;//切入点，切出点
while(true){//寻找切入边
esuc=suc[ein];
if(mv%(p[esuc]-p[ein])<0){
pin=Intersect(p[ein],p[esuc]-p[ein],m,mv);//传参顺序可能影响后面点积判断误差
if(pin==p[esuc]||(pin-p[ein])*(pin-p[esuc])/Dis(pin,p[esuc])<-EPS)
break;//点积判断要避免小数误差放大判断失效，加了个/Dis(pin,p[esuc])
}
ein=esuc;
if(ein==e)//未切入
return -1;
}
eout=suc[ein];
while(true){//寻找切出边
esuc=suc[eout];
if(mv%(p[esuc]-p[eout])>0){
pout=Intersect(p[eout],p[esuc]-p[eout],m,mv);
if(pout==p[esuc]||(pout-p[eout])*(pout-p[esuc])/Dis(pout,p[esuc])<-EPS)
break;//点积判断同上
}
eout=esuc;
if(eout==ein)//切于一个点，不算切入
return -1;
}
AddEdge(pin,pout,face[e],f);
p[eout]=pout;
suc[ein]=ecnt^1;
suc[ecnt^1]=pout==p[esuc]?esuc:eout;//处理退化情况
suc[ecnt]=ecnt-2;
ehead[face[e]]=ein;
return eout;
}
#define DEBUG(a,l,r) for(int i=l;i<=r;++i)cout<<#a<<'['<<i<<"]="<<a[i]<<endl;cout<<endl;
int main(){
freopen("in.in","r",stdin);
freopen("incremental.out","w",stdout);
int n,w,h;
cin>>n>>w>>h;

//初始化INF边界
a[n+1]=Point( INF, INF);
a[n+2]=Point(-INF, INF);
a[n+3]=Point(-INF,-INF);
a[n+4]=Point( INF,-INF);
AddEdge(Point(0,2*INF), Point(0,0),n+1,n+2);
AddEdge(Point(-2*INF,0),Point(0,0),n+2,n+3);
AddEdge(Point(0,-2*INF),Point(0,0),n+3,n+4);
AddEdge(Point(2*INF,0), Point(0,0),n+4,n+1);
ehead[n+1]=2,suc[2]=9,suc[9]=8;
ehead[n+2]=4,suc[4]=3,suc[3]=2;
ehead[n+3]=6,suc[6]=5,suc[5]=4;
ehead[n+4]=8,suc[8]=7,suc[7]=6;
//算法主体
for(int i=1;i<=n;++i){
cin>>a[i];
//寻找点i的最近点
int near;
double neardis=INF;
if(i==1)
near=n+(a[1].y>0?(a[1].x>0?1:2):(a[1].x<0?3:4));
else
for(int j=1;j<i;++j)
if(neardis>Dis(a[i],a[j]))
near=j,neardis=Dis(a[i],a[j]);
//构建点i的多边形
ehead[i]=ecnt+2;
int e=ehead[near];
do{
e=Cut((a[i]+a[face[e]])/2,(a[i]-a[face[e]]).Turn90(),e,i)^1;
}while(face[e]!=near);
suc[ehead[i]]=ecnt;
//		DEBUG(ehead,1,n+4);
//		DEBUG(face,2,ecnt);
//		DEBUG(suc,2,ecnt);
//		DEBUG(p,2,ecnt);
}
//加入长方形边界
double ans=0;
for(int i=1;i<=n;++i){
Cut(Point(0,0),Point(1,0),ehead[i],0);
Cut(Point(w,0),Point(0,1),ehead[i],0);
Cut(Point(w,h),Point(-1,0),ehead[i],0);
Cut(Point(0,h),Point(0,-1),ehead[i],0);
int cnt=0,e=ehead[i];
do{
if(ans<Dis(a[i],p[e]))ans=Dis(a[i],p[e]);
++cnt;
e=suc[e];
}while(e!=ehead[i]);
#ifdef TEST//测试生成多边形的边数，检查是否成功处理退化情况
cout<<cnt<<endl;
#endif
}
cout<<ans<<endl;
return 0;
}


## 点均匀随机分布前提下对逐点插入法的优化

//n=10000,time=280ms
#include<cstdio>
#include<cmath>
#include<iostream>
#include<algorithm>
#include<list>
#define double long double
#define TEST
using namespace std;
const int N=10009,E=50*N;
const double PI=acos(-1),EPS=1e-9,INF=1e4*(1+1/sqrt(2));
struct Point{
double x,y;
Point(){x=y=0;}
Point(double a){x=a;y=0;}
Point(double a,double b){x=a;y=b;}
friend istream&operator>>(istream&is,Point&a){return is>>a.x>>a.y;}
friend ostream&operator<<(ostream&os,Point a){return os<<a.x<<','<<a.y;}
Point operator-(){return Point(-x,-y);}
Point operator+(Point a){return Point(x+a.x,y+a.y);}
Point operator-(Point a){return Point(x-a.x,y-a.y);}
Point operator*(double a){return Point(x*a,y*a);}
Point operator/(double a){return Point(x/a,y/a);}
friend Point operator*(double a,Point b){b.x*=a,b.y*=a;return b;}
Point&operator+=(Point a){x+=a.x,y+=a.y;return*this;}
Point&operator-=(Point a){x-=a.x,y-=a.y;return*this;}
Point&operator*=(double a){x*=a,y*=a;return*this;}
Point&operator/=(double a){x/=a,y/=a;return*this;}
double operator*(Point a){return x*a.x+y*a.y;}
double operator%(Point a){return x*a.y-y*a.x;}
double Len(){return sqrt(x*x+y*y);}
friend double Dis(Point a,Point b){return (a-b).Len();}
Point Turn(double r){double c=cos(r),s=sin(r);return Point(x*c-y*s,y*c+x*s);}
Point Turn90(){return Point(-y,x);}
bool operator==(const Point&a){return fabs(x-a.x)<EPS&&fabs(y-a.y)<EPS;}
bool operator!=(const Point&a){return fabs(x-a.x)>EPS||fabs(y-a.y)>EPS;}
};
int ecnt=1,ehead[N],face[E],suc[E],id[N];
Point a[N],p[E];
double maxx[N];
list<int>li;

Point Intersect(Point a,Point av,Point b,Point bv){//直线求交点
return a+av*(bv%(b-a))/(bv%av);
}
void AddEdge(Point a,Point b,int af,int bf){//连边
p[++ecnt]=a;face[ecnt]=af;
p[++ecnt]=b;face[ecnt]=bf;
}
int Cut(Point m,Point mv,int e,int f){//有向直线(m,mv)切割边e所在面，另一侧划分给面f，返回翻边
int ein=e,eout,esuc;//切入边，切出边，后继
Point pin,pout;//切入点，切出点
while(true){//寻找切入边
esuc=suc[ein];
if(mv%(p[esuc]-p[ein])<0){
pin=Intersect(p[ein],p[esuc]-p[ein],m,mv);//传参顺序可能影响后面点积判断误差
if(pin==p[esuc]||(pin-p[ein])*(pin-p[esuc])/Dis(pin,p[esuc])<-EPS)
break;//点积判断要避免小数误差放大判断失效，加了个/Dis(pin,p[esuc])
}
ein=esuc;
if(ein==e)//未切入
return -1;
}
eout=suc[ein];
while(true){//寻找切出边
esuc=suc[eout];
if(mv%(p[esuc]-p[eout])>0){
pout=Intersect(p[eout],p[esuc]-p[eout],m,mv);
if(pout==p[esuc]||(pout-p[eout])*(pout-p[esuc])/Dis(pout,p[esuc])<-EPS)
break;//点积判断同上
}
eout=esuc;
if(eout==ein)//切于一个点，不算切入
return -1;
}
AddEdge(pin,pout,face[e],f);
p[eout]=pout;
suc[ein]=ecnt^1;
suc[ecnt^1]=pout==p[esuc]?esuc:eout;//处理退化情况
suc[ecnt]=ecnt-2;
ehead[face[e]]=ein;
maxx[face[e]]=0;//更新面最右侧x
e=ein;
do{
maxx[face[e]]=max(maxx[face[e]],p[e].x);
e=suc[e];
}while(e!=ein);
maxx[f]=max(maxx[f],pout.x);
return eout;
}
#define DEBUG(a,l,r) for(int i=l;i<=r;++i)cout<<#a<<'['<<i<<"]="<<a[i]<<endl;cout<<endl;
int main(){
freopen("in.in","r",stdin);
freopen("incremental1.out","w",stdout);
int n,w,h;
cin>>n>>w>>h;

//初始化INF边界
a[n+1]=Point( INF, INF);
a[n+2]=Point(-INF, INF);
a[n+3]=Point(-INF,-INF);
a[n+4]=Point( INF,-INF);
AddEdge(Point(0,2*INF), Point(0,0),n+1,n+2);
AddEdge(Point(-2*INF,0),Point(0,0),n+2,n+3);
AddEdge(Point(0,-2*INF),Point(0,0),n+3,n+4);
AddEdge(Point(2*INF,0), Point(0,0),n+4,n+1);
ehead[n+1]=2,suc[2]=9,suc[9]=8;
ehead[n+2]=4,suc[4]=3,suc[3]=2;
ehead[n+3]=6,suc[6]=5,suc[5]=4;
ehead[n+4]=8,suc[8]=7,suc[7]=6;
//算法主体
for(int i=1;i<=n;++i){
cin>>a[i];
id[i]=i;
}
sort(id+1,id+n+1,[](int x,int y){return a[x].x<a[y].x;});
for(int i=1;i<=n;++i){
//寻找点i的最近点
int near;
double neardis=INF;
if(i==1)
near=n+(a[1].y>0?(a[1].x>0?1:2):(a[1].x<0?3:4));
else{
auto iter=li.begin();
while(iter!=li.end()){
if(a[id[i]].x>maxx[*iter])
iter=li.erase(iter);//链表中去掉无用面，减少枚举
else{
if(neardis>Dis(a[id[i]],a[*iter]))
near=*iter,neardis=Dis(a[id[i]],a[*iter]);
++iter;
}
}
}
//构建点i的多边形
ehead[id[i]]=ecnt+2;
int e=ehead[near];
do{
e=Cut((a[id[i]]+a[face[e]])/2,(a[id[i]]-a[face[e]]).Turn90(),e,id[i])^1;
}while(face[e]!=near);
suc[ehead[id[i]]]=ecnt;
li.push_front(id[i]);
}
//加入长方形边界
double ans=0;
for(int i=1;i<=n;++i){
Cut(Point(0,0),Point(1,0),ehead[i],0);
Cut(Point(w,0),Point(0,1),ehead[i],0);
Cut(Point(w,h),Point(-1,0),ehead[i],0);
Cut(Point(0,h),Point(0,-1),ehead[i],0);
int cnt=0,e=ehead[i];
do{
if(ans<Dis(a[i],p[e]))ans=Dis(a[i],p[e]);
++cnt;
e=suc[e];
}while(e!=ehead[i]);
#ifdef TEST//测试生成多边形的边数，检查是否成功处理退化情况
cout<<cnt<<endl;
#endif
}
cout<<ans<<endl;
return 0;
}


# 分治法

函数merge(有序点集S)
如果|S|<=3
直接构造Voronoi图并返回
merge(S1)
merge(S2)
找到轮廓线的第一条边，p1p2的垂直平分线
遍历轮廓线，直到最后一条边
分别在S1,S2中做射线，与面p1、面p2产生交点
如果S1的交点更近
p1更新为交点所在边另一侧面对应基点
否则
p2更新为交点所在边另一侧面对应基点
在Voronoi图中连边
返回新Voronoi图


Samuel Peterson COMPUTING CONSTRAINED DELAUNAY TRIANGULATIONS

zball Delaunay剖分与平面欧几里得距离最小生成树

shunan 构造voronoi图分治算法

# 扫描线法

Seiyagoo Voronoi Diagram——维诺图

# 相关题目

LOJ6409「ICPC World Finals 2018」熊猫保护区

1 2 2
1 1 1

1.414

# 更新中！

posted @ 2021-09-17 01:53  Flash_Hu  阅读(430)  评论(4编辑  收藏  举报