# [ZJOI2018] bodyguard (Voronoi diagram) (Delaunay triangulation) (circle inversion) (3D convex hull) (Euler theorem)

Posted by wardo on Sun, 19 Jan 2020 05:56:19 +0100

# Portal

Please shout three times when you see this question: Du Jiao NB! Du Jiao NB! Du Jiao NB! In fact, it's almost the same whether Mr. Ji changes the data or not...

Anyway, no one wrote about the middle part

If Mr. Ji doesn't change the data range, it's estimated that the three-dimensional convex hull method of single logloglog can be found on the Internet now. I still don't single log three-dimensional convex hull

# Explanation:

For a point set of nnn size, the sum of the number of circumferences and triangles of any triangulation is 2n − 22n-22n − 2. The proof considers mapping it into a three-dimensional graph, and then proves it with the Euler theorem of topology.

Or directly consider using the Euler's theorem of plane graph (although it is essentially a thing), we have v − E+F=2V-E+F=2V − E+F=2

Let the number of points on the convex hull be kkk, it is obvious that the number of triangles is the number of faces, that is, F − 1F-1F − 1, the number of points is nnn, each edge on the convex hull is calculated twice, and each edge on the triangle is calculated three times, then

n−((F−1)∗3+k)/2+F=2n-((F-1)*3+k)/2+F=2n−((F−1)∗3+k)/2+F=2

Simplification is k+F − 1=2n − 2k+F-1=2n-2k+F − 1=2n − 2.

For convex hull, it is obvious that the triangulation required is convex.

So we can consider how many triangles there are in the case of convex periphery of triangulation.

It is easy to notice that the condition given by the problem is the form of inversion, which enlightens us to use the knowledge related to inversion to solve it.

The circle is the most closely related to inversion, and the triangulation which is the most closely related to circle is section D.

Define the inner circle as a circle passing through at least three points in the point set without any other points inside, and the outer circle as a circle passing through at least three points in the point set without any other points outside.

It is easy to find that the inner circle can be directly calculated by the nearest V-graph, that is, directly up the D-section. The three points of the outer circle correspond to the three adjacent points in the farthest point V graph. At present, it is not clear whether there is a fast direct method

Considering the graph relationship before and after inversion, it is obvious that:

1. For the inner circle, if the inversion center is in the inner circle, its inversion is the outer circle, otherwise it is still the inner circle.
2. For the outer circle, if the inversion center is in the inner circle, its inversion is the inner circle, otherwise it is still the outer circle.

The probability of inversion center on the boundary is 000, not considered.

So we can find out all the inner and outer circles, and then find the intersection with the original rectangle.

As the outer circle can't be found, it seems to be in a stalemate.

Take a look at the equation of a circle: (x − x0)2+(y − y0) 2=R2 (x-x ^ 0) ^ 2+(y-Y ^ 0) ^ 2 = R ^ 2 (x − x0)2+(y − y0) 2=R2

The form of splitting is as follows: x2+y2+Ax+By+C=0x^2+y^2+Ax+By+C=0x2+y2+Ax+By+C=0.

The plane equation in three-dimensional coordinate system is z+Ax+By+C=0z+Ax+By+C=0z+Ax+By+C=0.

Consider mapping the points (x,y)(x,y)(x,y) of the two-dimensional plane to the three-dimensional (x,y,x2+y2)(x,y,x^2+y^2)(x,y,x2+y2), so that the plane of each three points can represent a circle. Next, consider how to represent the relationship between a circle and a point. Obviously, if the point is below the plane, it is inside the circle, otherwise it is outside the circle.

So we can directly find the three-dimensional convex hull, which faces the negative direction of zzz axis to represent the inner circle, and faces the positive direction to represent the outer circle.

In order to avoid parallel to zzz axis, a little random disturbance can be added.

Code:

```#include<bits/stdc++.h>
#define ll long long
#define db long double
#define re register
#define cs const

namespace IO{
inline char gc(){
static cs int Rlen=1<<22|1;static char buf[Rlen],*p1,*p2;
}template<typename T>T get_integer(){
char c;T num;while(!isdigit(c=gc()));num=c^48;
while(isdigit(c=gc()))num=((num+(num<<2))<<1)+(c^48);return num;
}inline int gi(){return get_integer<int>();}
}using namespace IO;

using std::cerr;
using std::cout;

std::mt19937 R(time(0));
std::uniform_int_distribution<int> rnd(0,(int)1e9);
cs db eps=5e-8;
inline db rndb(){return rnd(R)/1e10;}
inline db reps(){return (rndb()-.5)*eps;}

inline int sign(db x){return x<-eps?-1:(x>eps?1:0);}
struct Pnt{
db x,y;Pnt(){}Pnt(db _x,db _y):x(_x),y(_y){}
friend Pnt operator+(cs Pnt &a,cs Pnt &b){return Pnt(a.x+b.x,a.y+b.y);}
friend Pnt operator-(cs Pnt &a,cs Pnt &b){return Pnt(a.x-b.x,a.y-b.y);}
friend db operator*(cs Pnt &a,cs Pnt &b){return a.x*b.y-b.x*a.y;}
friend db dot(cs Pnt &a,cs Pnt &b){return a.x*b.x+a.y*b.y;}
friend Pnt operator*(cs Pnt &a,db b){return Pnt(a.x*b,a.y*b);}
db len()cs{return sqrtl(x*x+y*y);}db ang()cs{return atan2l(y,x);}
};

struct pt{
db x,y,z;pt(){}pt(db _x,db _y,db _z):x(_x),y(_y),z(_z){}
friend pt operator+(cs pt &a,cs pt &b){return pt(a.x+b.x,a.y+b.y,a.z+b.z);}
friend pt operator-(cs pt &a,cs pt &b){return pt(a.x-b.x,a.y-b.y,a.z-b.z);}
friend pt operator*(cs pt &a,cs pt &b){return pt(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x);}
friend db dot(cs pt &a,cs pt &b){return a.x*b.x+a.y*b.y+a.z*b.z;}
friend bool operator==(cs pt &a,cs pt &b){return !sign(a.x-b.x)&&!sign(a.y-b.y)&&!sign(a.z-b.z);}
friend bool operator!=(cs pt &a,cs pt &b){return sign(a.x-b.x)||sign(a.y-b.y)||sign(a.z-b.z);}
db len()cs{return sqrtl(x*x+y*y+z*z);}
};

inline db crs(cs Pnt &a,cs Pnt &b,cs Pnt &c){
return (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x);
}
inline pt crs(cs pt &a,cs pt &b,cs pt &c){
return (b-a)*(c-a);
}
inline pt trans(cs Pnt &a){
return pt(a.x,a.y,a.x*a.x+a.y*a.y);
}
inline Pnt trans(cs pt &a){
return Pnt(a.x,a.y);
}

inline db area(cs pt &a,cs pt &b,cs pt &c){
return crs(a,b,c).len();
}
inline db vol(cs pt &a,cs pt &b,cs pt &c,cs pt &d){
return dot(crs(a,b,c),d-a);
}

cs int N=2e3+7;

int n,ct,id[N][N];
bool ok[N<<2|1];

pt p[N];Pnt q[N];
db xa,xb,ya,yb,S;
Pnt poly,tmp;

struct Face{
int a,b,c;Face(){}Face(int _a,int _b,int _c):a(_a),b(_b),c(_c){}
pt normal()cs{return crs(p[a],p[b],p[c]);}
db area()cs{return normal().len();}
db vol(cs pt &d)cs{return dot(normal(),d-p[a]);}
}f[N<<2|1];

void dfs(int,int);
void deal(int,int,int);

void deal(int x,int a,int b){
if(ok[id[a][b]]){
if(f[id[a][b]].vol(p[x])>0)dfs(x,id[a][b]);
else {
id[b][a]=id[a][x]=id[x][b]=++ct;
f[ct]=Face(b,a,x);ok[ct]=true;
}
}
}

void dfs(int p,int x){
ok[x]=false;
deal(p,f[x].b,f[x].a);
deal(p,f[x].a,f[x].c);
deal(p,f[x].c,f[x].b);
}

void deal(int p,int x){
ok[x]=false;
deal(p,f[x].b,f[x].a);
deal(p,f[x].a,f[x].c);
deal(p,f[x].c,f[x].b);
}

void Convex_3D(){
ct=0;int flag=true;
for(int re i=2;i<=n;++i)if(p!=p[i])
{std::swap(p,p[i]);flag=false;break;}
if(flag)printf("%.20Lf",(db)1),exit(0);flag=true;
for(int re i=3;i<=n;++i)if(sign(area(p,p,p[i])))
{std::swap(p,p[i]);flag=false;break;}
if(flag)printf("%.20Lf",(db)n),exit(0);flag=true;
for(int re i=4;i<=n;++i)if(sign(vol(p,p,p,p[i])))
{std::swap(p,p[i]);flag=false;break;}
if(flag)printf("%.20Lf",(db)n),exit(0);

for(int re i=1;i<=n;++i)
p[i].x+=reps(),p[i].y+=reps(),p[i].z+=reps();
for(int re i=1;i<=4;++i){
int a=(i&3)+1,b=(a&3)+1,c=(b&3)+1;
if(vol(p[a],p[b],p[c],p[i])>0)std::swap(b,c);
id[a][b]=id[b][c]=id[c][a]=++ct;
f[ct]=Face(a,b,c);ok[ct]=true;
}
for(int re i=5;i<=n;++i){
for(int re j=1;j<=ct;++j)
if(f[j].vol(p[i])>0)
{dfs(i,j);flag=j;break;}
if(!flag)continue;int cur=flag-1;flag=false;
for(int re j=cur+1;j<=ct;++j)if(ok[j]){
if((++cur)==j)continue;f[cur]=f[j],ok[cur]=true;
id[f[cur].a][f[cur].b]=id[f[cur].b][f[cur].c]=id[f[cur].c][f[cur].a]=cur;
}ct=cur;
}
}

inline db solve(db r,db x,db y){
if(x<0||y<0)return (x*y<0)?-solve(r,fabsl(x),fabsl(y)):solve(r,fabsl(x),fabsl(y));
if(x*x+y*y<=r*r)return x*y;db p=(r<y)?0:sqrtl(r*r-y*y),q=(r<x)?0:sqrtl(r*r-x*x);
return (p*y+q*x+(atan2l(y,p)-atan2(q,x))*r*r)*0.5;
}
inline Pnt intersection(cs Pnt &a,cs Pnt &b,cs Pnt &c,cs Pnt &d){
return a+(b-a)*(crs(c,d,a)/((b-a)*(d-c)));
}
inline bool on_seg(cs Pnt &p,cs Pnt &s,cs Pnt &e){return sign(dot(p-s,p-e))<=0;}

db calc(cs Face &f){
int a=f.a,b=f.b,c=f.c;std::swap(b,c);
if(fabsl(crs(q[a],q[b],q[c]))<eps){
db res=0;int tot=0;
if((crs(p[a],p[b],p[c]).z<0)^
(dot(p[b]-p[a],p[c]-p[a])<0))
std::swap(b,c);
for(int re i=0;i<4;++i){
if(crs(q[b],q[c],poly[i])>=-eps)tmp[tot++]=poly[i];
if(fabsl((q[c]-q[b])*(poly[i]-poly[i+1]))>eps){
Pnt p=intersection(q[b],q[c],poly[i],poly[i+1]);
if(on_seg(p,poly[i],poly[i+1]))tmp[tot++]=p;
}
}for(int re i=1;i+1<tot;++i)res+=crs(tmp,tmp[i],tmp[i+1]);
return fabsl(res*0.5);
}pt P=crs(p[a],p[b],p[c]);
db x=P.x/P.z*-.5,y=P.y/P.z*-.5,r=sqrtl((p[a].x-x)*(p[a].x-x)+(p[a].y-y)*(p[a].y-y));
db xl=std::max(xa-x,-r),xr=std::min(xb-x,r),yd=std::max(ya-y,-r),yu=std::min(yb-y,r);
if(xl>=xr||yd>=yu)return 0;
return solve(r,-xl,-yd)+solve(r,xr,-yd)+solve(r,-xl,yu)+solve(r,xr,yu);
}

void Main(){
n=gi();if(n==3)return (void)printf("3\n");
xa=gi(),ya=gi(),xb=gi(),yb=gi();
poly=Pnt(xa,ya),poly=Pnt(xb,ya);
poly=Pnt(xb,yb),poly=Pnt(xa,yb);poly=poly;
S=(xb-xa)*(yb-ya);db s=(2*n-2)*S;
for(int re i=1;i<=n;++i)
q[i].x=gi(),q[i].y=gi(),p[i]=trans(q[i]);
Convex_3D();
for(int re i=1;i<=ct;++i)
s-=(f[i].normal().z>0)?calc(f[i]):(S-calc(f[i]));
printf("%.20Lf",s/S);
}

void file(){
#ifdef zxyoi
freopen("guard.in","r",stdin);
#endif
}
signed main(){file();Main();return 0;}
```  967 original articles published, 369 praised, 70000 visitors+