# 半平面交

#### ( 一 )求解半平面交

Uyuw's Concert

``````#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
const int MAXN=40010;
const double EPS=1e-10;
struct Point
{
double x,y;
Point(double x=0,double y=0):x(x),y(y){}
};
typedef Point Vector;
Vector operator -(Vector A,Vector B)
{
return Vector(A.x-B.x,A.y-B.y);
}
Vector operator+(Vector A,Vector B)
{
return Vector(A.x+B.x,A.y+B.y);
}
double cross(Vector A,Vector B)
{
return A.x*B.y-A.y*B.x;
}
Vector operator*(Vector A,double val)
{
return Vector(A.x*val,A.y*val);
}
struct Line
{
Point p;//直线上的任意一点
Vector v;//方向向量，它的左边是对应的半平面
double ang;//极角
Line(){}
Line(Point p,Vector v):p(p),v(v){  ang=atan2(v.y,v.x); }
bool operator<(const Line &l) const
{
return ang<l.ang;
}
};
bool onLeft(Line L,Point p)//点p在有向直线的左边(在线上不算)
{
return cross(L.v,p-L.p)>0;
}
Point getIntersection(Line a,Line b)//两直线交点，假定交点唯一存在
{
Vector u=a.p-b.p;
double t=cross(b.v,u)/cross(a.v,b.v);
return a.p+a.v*t;
}
Point p[MAXN];
Line q[MAXN];//双端队列
Line line[MAXN];
Point out[MAXN];
//半平面交的过程
int halfInter(Line *L,int n,Point *poly)
{
sort(L,L+n);
int first,last;
q[first=last=0]=L[0];
for(int i=1;i<n;i++)
{
while(first<last&&!onLeft(L[i],p[last-1])) last--;
while(first<last&&!onLeft(L[i],p[first])) first++;
q[++last]=L[i];
if(fabs(cross(q[last].v,q[last-1].v))<EPS)
{//两直线平行,取内侧的一个
last--;
if(onLeft(q[last],L[i].p))
{
q[last]=L[i];
}
}
if(first<last) p[last-1]=getIntersection(q[last],q[last-1]);
}
while(first<last&&!onLeft(q[first],p[last-1])) last--;//删除无用平面
if(last-first<=1) return 0;//空集
p[last]=getIntersection(q[last],q[first]);//计算首尾两个半平面的交点
int m=0;
for(int i=first;i<=last;i++)
{
poly[m++]=p[i];
}
return m;
}
double polyArea(Point *p,int n)
{
double area=0.0;
for(int i=1;i<n-1;i++)
{
area+=cross(p[i]-p[0],p[i+1]-p[0]);
}
return area/2;
}
int main()
{
int n,m=0;
scanf("%d",&n);
double x1,y1,x2,y2;
for(int i=0;i<n;i++)
{
scanf("%lf%lf%lf%lf",&x1,&y1,&x2,&y2);
line[m++]=Line(Point(x1,y1),Vector(x2-x1,y2-y1));
}
for(int i=0;i<4;i++)
{
line[m++]=Line(Point(0,0),Vector(1,0));
line[m++]=Line(Point(10000,0),Vector(0,1));
line[m++]=Line(Point(10000,10000),Vector(-1,0));
line[m++]=Line(Point(0,10000),Vector(0,-1));
}
int res=halfInter(line,m,out);
double area=polyArea(out,res);
printf("%.1f\n",area);
return 0;
}
``````

hdu 1632 Polygons

#### ( 二 )求解多边形的核

Rotating Scoreboard

``````#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
const int MAXN=1010;
const double EPS=1e-10;
struct Point
{
double x,y;
Point(double x=0,double y=0):x(x),y(y){}
};
typedef Point Vector;
Vector operator-(Vector A,Vector B)
{
return Vector(A.x-B.x,A.y-B.y);
}
Vector operator+(Vector A,Vector B)
{
return Vector(A.x+B.x,A.y+B.y);
}
Vector operator *(Vector A,double val)
{
return Vector(A.x*val,A.y*val);
}
double cross(Vector A,Vector B)
{
return A.x*B.y-A.y*B.x;
}
struct Line
{
Point p;
Vector v;
double ang;
Line(){}
Line(Point p,Vector v):p(p),v(v){ ang=atan2(v.y,v.x); }
bool operator<(const Line &l) const
{
return ang<l.ang;
}
};
int dcmp(double val)
{
if(abs(val)<EPS) return 0;
return val<0?-1:1;
}
bool onLeft(Line line,Point p)//点在直线上也可以,同时考虑下误差
{
return dcmp(cross(line.v,p-line.p))>=0;
}
Point getIntersection(Line a,Line b)
{
Vector u=a.p-b.p;
double t=cross(b.v,u)/cross(a.v,b.v);
return a.p+a.v*t;
}
Line line[MAXN];
Line que[MAXN];
Point p[MAXN];
Point ans[MAXN];
Point point[MAXN];
int halfIntersection(Line *L,int n,Point *poly)
{
sort(L,L+n);
int first,last;
que[first=last=0]=L[0];
for(int i=1;i<n;i++)
{
while(first<last&&!onLeft(L[i],p[last-1])) last--;
while(first<last&&!onLeft(L[i],p[first])) first++;
que[++last]=L[i];
if(dcmp(cross(que[last].v,que[last-1].v))==0)
{
last--;
if(onLeft(que[last],L[i].p)) que[last]=L[i];
}
if(first<last) p[last-1]=getIntersection(que[last],que[last-1]);
}
while(first<last&&!onLeft(que[first],p[last-1])) last--;
if(last-first<=1) return 0;
p[last]=getIntersection(que[first],que[last]);
int m=0;
for(int i=first;i<=last;i++)
{
poly[m++]=p[i];
}
return m;
}
int main()
{
int n,t,res;
double x,y;
scanf("%d",&t);
while(t--)
{
scanf("%d",&n);
for(int i=0;i<n;i++)
{
scanf("%lf%lf",&point[i].x,&point[i].y);
}
point[n]=point[0];
int cnt=0;
for(int i=1;i<=n;i++)
{
line[cnt++]=Line(point[i-1],point[i-1]-point[i]);
}
res=halfIntersection(line,cnt,ans);
if(res) printf("YES\n");
else printf("NO\n");
}
return 0;
}
``````

#### ( 三 )求解凸多边形的最大内切圆

Most Distant Point from the Sea

``````#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
const int MAXN=1510;
const double EPS=1e-10;
struct Point
{
double x,y;
Point(double x=0,double y=0):x(x),y(y){}
};
typedef Point Vector;
Vector operator -(Vector A,Vector B)
{
return Vector(A.x-B.x,A.y-B.y);
}
Vector operator +(Vector A,Vector B)
{
return Vector(A.x+B.x,A.y+B.y);
}
Vector operator*(Vector A,double val)
{
return Vector(A.x*val,A.y*val);
}
double cross(Vector A,Vector B)
{
return A.x*B.y-A.y*B.x;
}
struct Line
{
Point p;
Vector v;
double ang;
Line(){}
Line(Point p,Vector v):p(p),v(v){ ang=atan2(v.y,v.x); }
bool operator <(const Line &l)const
{
return ang<l.ang;
}
};
int dcmp(double val)
{
if(abs(val)<EPS) return 0;
return val>0?1:-1;
}
bool onLeft(Line line,Point p)
{
return dcmp(cross(line.v,p-line.p))>=0;
}
Point getIntersection(Line a,Line b)
{
Vector u=a.p-b.p;
double t=cross(b.v,u)/cross(a.v,b.v);
return a.p+a.v*t;
}
Point p[MAXN];
Point point[MAXN];
Point change[MAXN];
Line que[MAXN];
Line line[MAXN];
Vector normal[MAXN];
bool halfIntersection(Line *L,int n)
{
sort(L,L+n);
int first,last;
que[first=last=0]=L[0];
for(int i=1;i<n;i++)
{
while(first<last&&!onLeft(L[i],p[last-1])) last--;
while(first<last&&!onLeft(L[i],p[first])) first++;
que[++last]=L[i];
if(dcmp(cross(que[last].v,que[last-1].v))==0)
{
last--;
if(onLeft(que[last],L[i].p)) que[last]=L[i];
}
if(first<last) p[last-1]=getIntersection(que[last],que[last-1]);
}
while(first<last&&!onLeft(que[first],p[last-1])) last--;
if(last-first<=1) return false;
else return true;
}
double polyArea(Point *p,int n)
{
double area=0;
for(int i=1;i<n-1;i++)
{
area+=cross(p[i]-p[0],p[i+1]-p[0]);
}
return area/2;
}
double length(Vector A)
{
return sqrt(A.x*A.x+A.y*A.y);
}
Vector Normal(Vector A)//单位法向量
{
double len=length(A);
return Vector(-A.y/len,A.x/len);
}
int main()
{
int n;
double x,y,lef,rig,mid;
while(scanf("%d",&n)!=EOF,n)
{
for(int i=0;i<n;i++)
{
scanf("%lf%lf",&point[i].x,&point[i].y);
}
point[n]=point[0];
for(int i=1;i<=n;i++)
{
normal[i-1]=Normal(point[i]-point[i-1]);
}
lef=0.0;rig=10000.0;
while(lef+EPS<rig)
{
mid=(lef+rig)/2.0;
for(int i=1;i<=n;i++)
{
line[i-1]=Line(point[i-1]+normal[i-1]*mid,point[i]-point[i-1]);//直线向垂直方向平移
}
if(halfIntersection(line,n)) lef=mid;
else rig=mid;
}
printf("%.6f\n",mid);
}
return 0;
}
``````

Feng Shui

``````#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
const int MAXN=1510;
const double EPS=1e-10;
int dcmp(double val)
{
if(abs(val)<EPS) return 0;
return val>0?1:-1;
}
struct Point
{
double x,y;
Point(double x=0,double y=0):x(x),y(y){}
};
typedef Point Vector;
Vector operator -(Vector A,Vector B)
{
return Vector(A.x-B.x,A.y-B.y);
}
Vector operator +(Vector A,Vector B)
{
return Vector(A.x+B.x,A.y+B.y);
}
Vector operator*(Vector A,double val)
{
return Vector(A.x*val,A.y*val);
}
double cross(Vector A,Vector B)
{
return A.x*B.y-A.y*B.x;
}
struct Line
{
Point p;
Vector v;
double ang;
Line(){}
Line(Point p,Vector v):p(p),v(v){ ang=atan2(v.y,v.x); }
bool operator <(const Line &l)const
{
return ang<l.ang;
}
};
bool onLeft(Line line,Point p)
{
return dcmp(cross(line.v,p-line.p))>=0;
}
Point getIntersection(Line a,Line b)
{
Vector u=a.p-b.p;
double t=cross(b.v,u)/cross(a.v,b.v);
return a.p+a.v*t;
}
Point p1,p2;
Point p[MAXN];
Point point[MAXN];
Point ans[MAXN];
Line que[MAXN];
Line line[MAXN];
int halfIntersection(Line *L,int n,Point *poly)
{
sort(L,L+n);
int first,last;
que[first=last=0]=L[0];
for(int i=1;i<n;i++)
{
while(first<last&&!onLeft(L[i],p[last-1])) last--;
while(first<last&&!onLeft(L[i],p[first])) first++;
que[++last]=L[i];
if(dcmp(cross(que[last].v,que[last-1].v))==0)
{
last--;
if(onLeft(que[last],L[i].p)) que[last]=L[i];
}
if(first<last) p[last-1]=getIntersection(que[last],que[last-1]);
}
while(first<last&&!onLeft(que[first],p[last-1])) last--;
if(last-first<=1) return 0;
p[last]=getIntersection(que[first],que[last]);
int m=0;
for(int i=first;i<=last;i++)
{
poly[m++]=p[i];
}
return m;
}
double length(Vector A)
{
return sqrt(A.x*A.x+A.y*A.y);
}
Vector normal(Vector A)
{
double len=length(A);
return Vector(-A.y/len,A.x/len);
}
void rotateCalipers(Point *ch,int n)
{
double res=-1,area,len;
int up=1;
ch[n]=ch[0];
for(int i=1;i<=n;i++)
{
while(cross(ch[i]-ch[i-1],ch[up+1]-ch[i-1])>cross(ch[i]-ch[i-1],ch[up]-ch[i-1])) up=(up+1)%n;
len=length(ch[up+1]-ch[i]);
if(dcmp(len-res)>=0)
{
p1=ch[up+1];
p2=ch[i];
res=len;
}
len=length(ch[up]-ch[i-1]);
if(dcmp(len-res)>=0)
{
p1=ch[i-1];
p2=ch[up];
res=len;
}
}
}
int main()
{
int n,mid;
double x,y,r,maxr,d;
while(scanf("%d%lf",&n,&r)!=EOF)
{
maxr=0.0;
for(int i=0;i<n;i++)
{
scanf("%lf%lf",&point[i].x,&point[i].y);
}
point[n]=point[0];
for(int i=1;i<=n;i++)
{
line[i-1]=Line(point[i-1]+normal(point[i-1]-point[i])*r,point[i-1]-point[i]);
}
int res=halfIntersection(line,n,ans);
rotateCalipers(ans,res);
printf("%lf %lf %lf %lf\n",p1.x,p1.y,p2.x,p2.y);
}
return 0;
}
``````