【kuangbin专题】计算几何_半平面交

时间:2023-03-08 21:07:54
【kuangbin专题】计算几何_半平面交

1.poj3335 Rotating Scoreboard

传送:http://poj.org/problem?id=3335

题意:就是有个球场,球场的形状是个凸多边形,然后观众是坐在多边形的边上的,问你是否在球场上有个地方可以放一个记分牌,然后所有的观众都可以看得到的。

分析:多边形是否存在内核。

 #include<iostream>
#include<algorithm>
#include<cstring>
#include<cmath>
using namespace std;
const int maxn=;
const double eps=1e-;
int sgn(double x){
if (fabs(x)<eps) return ;
if (x<) return -;
return ;
}
struct point{
double x,y;
point(){}
point(double _x,double _y):x(_x),y(_y){}
point operator -(const point &b)const{
return point(x-b.x,y-b.y);
}
double operator ^(const point &b)const{
return x*b.y-y*b.x;
}
}p[maxn];
struct line{
point s,e;
line(){}
line(point _s,point _e):s(_s),e(_e){}
point crosspoint(line v){
int a1=(v.e-v.s)^(s-v.s),a2=(v.e-v.s)^(e-v.s);
return point((s.x*a2-e.x*a1)/(a2-a1),(s.y*a2-e.y*a1)/(a2-a1));
}
//两向量平行(对应直线平行或重合)
bool parallel(line v){
return sgn((e-s)^(v.e-v.s))==;
}
};
struct halfplane:public line{
double angle;
halfplane(){}
halfplane(line v){s=v.s;e=v.e;}
void calcangle(){
angle=atan2(e.y-s.y,e.x-s.x);
}
bool operator <(const halfplane &b)const{
return angle<b.angle;
}
};
struct halfplanes{
int n;
halfplane hp[maxn];
point p[maxn];
int que[maxn]; int st,ed; //双端队列
void push(halfplane tmp){hp[n++]=tmp;}
//去重
void unique(){
int m=;
for (int i=;i<n;i++){
if (sgn(hp[i].angle-hp[i-].angle)!=) hp[m++]=hp[i];
else if (sgn((hp[m-].e-hp[m-].s)^(hp[i].s-hp[m-].s))>) hp[m-]=hp[i];
}
n=m;
}
bool halfplaneinsert(){
for (int i=;i<n;i++) hp[i].calcangle();
sort(hp,hp+n);
unique();
que[st=]=; que[ed=]=;
p[]=hp[].crosspoint(hp[]);
for (int i=;i<n;i++){
while (st<ed && sgn((hp[i].e-hp[i].s)^(p[ed]-hp[i].s))<) ed--;
while (st<ed && sgn((hp[i].e-hp[i].s)^(p[st+]-hp[i].s))<) st++;
que[++ed]=i;
if (hp[i].parallel(hp[que[ed-]])) return false;
p[ed]=hp[i].crosspoint(hp[que[ed-]]);
}
while (st<ed && sgn((hp[que[st]].e-hp[que[st]].s)^(p[ed]-hp[que[st]].s))<) ed--;
while (st<ed && sgn((hp[que[ed]].e-hp[que[ed]].s)^(p[st+]-hp[que[ed]].s))<) st++;
if (st+>=ed) return false; else return true;
}
};
halfplanes ha;
int main(){
int t,n; cin >>t;
while (t--){
cin >> n;
ha.n=;
for (int i=;i<n;i++) cin >> p[i].x >> p[i].y;
for (int i=;i<n;i++){
ha.push(halfplane(line(p[i],p[(i-+n)%n])));
}
if (ha.halfplaneinsert()) cout << "YES\n"; else cout << "NO\n";
}
return ;
}

poj3335

2.poj1279 Art Gallery

传送:http://poj.org/problem?id=1279

题意:给一个多边形(顺时针给出点),求内核面积。

分析:rt。

 #include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const int maxn=;
const double eps=1e-;
int sgn(double x){
if (fabs(x)<eps) return ;
if (x<) return -;
return ;
}
struct point{
double x,y;
point(){}
point(double _x,double _y):x(_x),y(_y){}
point operator -(const point &b)const{
return point(x-b.x,y-b.y);
}
double operator ^(const point &b)const{
return x*b.y-y*b.x;
}
}p[maxn];
struct line{
point s,e;
line(){}
line(point _s,point _e):s(_s),e(_e){}
point crosspoint(line v){
int a1=(v.e-v.s)^(s-v.s),a2=(v.e-v.s)^(e-v.s);
return point((s.x*a2-e.x*a1)/(a2-a1),(s.y*a2-e.y*a1)/(a2-a1));
}
//两向量平行(对应直线平行或重合)
bool parallel(line v){
return sgn((e-s)^(v.e-v.s))==;
}
};
struct polygon{
int n;
point p[maxn];
double getarea(){
double sum=;
for (int i=;i<n;i++) sum+=(p[i]^(p[(i+)%n]));
return fabs(sum)/;
}
};
struct halfplane:public line{
double angle;
halfplane(){}
halfplane(line v){s=v.s;e=v.e;}
void calcangle(){
angle=atan2(e.y-s.y,e.x-s.x);
}
bool operator <(const halfplane &b)const{
return angle<b.angle;
}
};
polygon C;
struct halfplanes{
int n;
halfplane hp[maxn];
point p[maxn];
int que[maxn]; int st,ed; //双端队列
void push(halfplane tmp){hp[n++]=tmp;}
//去重
void unique(){
int m=;
for (int i=;i<n;i++){
if (sgn(hp[i].angle-hp[i-].angle)!=) hp[m++]=hp[i];
else if (sgn((hp[m-].e-hp[m-].s)^(hp[i].s-hp[m-].s))>) hp[m-]=hp[i];
}
n=m;
}
bool halfplaneinsert(){
for (int i=;i<n;i++) hp[i].calcangle();
sort(hp,hp+n);
unique();
que[st=]=; que[ed=]=;
p[]=hp[].crosspoint(hp[]);
for (int i=;i<n;i++){
while (st<ed && sgn((hp[i].e-hp[i].s)^(p[ed]-hp[i].s))<) ed--;
while (st<ed && sgn((hp[i].e-hp[i].s)^(p[st+]-hp[i].s))<) st++;
que[++ed]=i;
if (hp[i].parallel(hp[que[ed-]])) return false;
p[ed]=hp[i].crosspoint(hp[que[ed-]]);
}
while (st<ed && sgn((hp[que[st]].e-hp[que[st]].s)^(p[ed]-hp[que[st]].s))<) ed--;
while (st<ed && sgn((hp[que[ed]].e-hp[que[ed]].s)^(p[st+]-hp[que[ed]].s))<) st++;
if (st+>=ed) return false; else return true;
}
//得到最后半平面交的凸多边形
//得到先调用halfplaneinsert(),且返回true
void getconvex(polygon &con){
p[st]=hp[que[st]].crosspoint(hp[que[ed]]);
con.n=ed-st+;
for (int j=st,i=;j<=ed;j++,i++) con.p[i]=p[j];
}
};
halfplanes ha;
int main(){
int t,n; cin >> t;
while (t--){
cin >> n;
for (int i=;i<n;i++) cin >> p[i].x >> p[i].y;
ha.n=;
for (int i=;i<n;i++) ha.push(line(p[i],p[(i-+n)%n]));
ha.halfplaneinsert();
C.n=;
ha.getconvex(C);
printf("%.2f\n",C.getarea());
}
return ;
}

poj1279

3.poj3525 Most Distant Point from the Sea

传送:http://poj.org/problem?id=3525

题意:给出一个海岛,求出岛上离海最远的点。输出距离。

分析:半平面交+二分。就是求多边形的最大半径圆   //板子

 #include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const double eps=1e-;
const int maxn=;
inline double sqr(double x){return x*x;}
int N;
int sgn(double x){
if (fabs(x)<eps) return ;
if (x<) return -;
return ;
}
struct point{
double x,y;
point(){}
point(double _x,double _y):x(_x),y(_y){}
point operator -(const point &b)const{
return point(x-b.x,y-b.y);
}
double operator ^(const point &b)const{
return x*b.y-y*b.x;
}
};
struct line{
point s,e;
double k; //斜率
line(){}
line(point _s,point _e):s(_s),e(_e){k=atan2(e.y-s.y,e.x-s.x);}
point operator &(const line &b)const{ //求两直线交点
point res=s;
double t=((s-b.s)^(b.s-b.e))/((s-e)^(b.s-b.e));
res.x+=(e.x-s.x)*t; res.y+=(e.y-s.y)*t;
return res;
}
};
point p[maxn]; //记录最初给的点集
line L[maxn]; //由最初的点集生成直线的集合
point pp[maxn]; //记录半平面交的结果的点集
//半平面交,直线的左边代表有效区域
bool HPIcmp(line a,line b){ //直线排序函数
if(sgn(a.k-b.k)!=) return a.k<b.k; //斜率排序
//斜率相同我也不知道怎么办
return ((a.s-b.s)^(b.e-b.s))<;
}
void HPI(line L[],int n,point res[],int &resn){
//line是半平面交的直线的集合 n是直线的条数 res是结果
//的点集 resn是点集里面点的个数
int tot=n;
sort(L,L+n,HPIcmp);
tot=;
//去掉斜率重复的
for(int i=;i<n;i++) if(sgn(L[i].k-L[i-].k)!=) L[tot++]=L[i];
int head,tail; line Q[maxn]; //双端队列
Q[head=]=L[]; Q[tail=]=L[];
resn=;
for(int i=;i<tot;i++){
if(sgn((Q[tail].e-Q[tail].s)^(Q[tail-].e-Q[tail-].s))==||
sgn((Q[head].e-Q[head].s)^(Q[head+].e-Q[head+].s))==) return;
while(head<tail && (((Q[tail]&Q[tail-])-L[i].s)^(L[i].e-L[i].s))>eps) tail--;
while(head<tail && (((Q[head]&Q[head+])-L[i].s)^(L[i].e-L[i].s))>eps) head++;
Q[++tail]=L[i];
}
while(head<tail && (((Q[tail]&Q[tail-])-Q[head].s)^(Q[head].e-Q[head].s))>eps) tail--;
while(head<tail && (((Q[head]&Q[head-])-Q[tail].s)^(Q[tail].e-Q[tail].e))>eps) head++;
if(tail<=head+) return;
for(int i=head;i<tail; i++) res[resn++]=Q[i]&Q[i+];
if(head<tail-) res[resn++]=Q[head]&Q[tail];
}
double dist(point a,point b){
return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y));
}
void change(point p1,point p2,point &a,point &b,double p){
double len=dist(p1,p2);
double dx=(p1.y-p2.y)*p/len,dy=(p2.x-p1.x)*p/len;
a.x=p1.x+dx;a.y=p1.y+dy;
b.x=p2.x+dx;b.y=p2.y+dy;
}
bool check(double mid){
for (int i=;i<N;i++){
point a,b;
change(p[i],p[(i+)%N],a,b,mid);
L[i]=line(a,b);
}
int resn;
HPI(L,N,pp,resn); //resn等于0说明移多了
if (resn==) return false; else return true;
}
int main(){
while (cin >> N && N){
for (int i=;i<N;i++)
cin >> p[i].x >> p[i].y;
double l=,r=,mid;
double ans;
while (r-l>=eps){
mid=(l+r)/;
if (check(mid)){
ans=mid;
l=mid+eps;
}
else r=mid-eps;
}
printf("%.6f\n",ans);
}
return ;
}

poj3525

4.poj3384 Feng Shui

传送:http://poj.org/problem?id=3384

题意:给定两个同样大小的圆放置在多边形内(可以与多边形相切,两圆可以重叠),尽可能大的覆盖多边形面积,问两圆圆心。

分析:用半平面交将多边形的每条边一起向“内”推进R,得到新的多边形,然后求多边形的最远两点。

 #include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const double eps=1e-;
const int maxn=;
inline double sqr(double x){return x*x;}
int n; double r;
int sgn(double x){
if (fabs(x)<eps) return ;
if (x<) return -;
return ;
}
struct point{
double x,y;
point(){}
point(double _x,double _y):x(_x),y(_y){}
point operator -(const point &b)const{
return point(x-b.x,y-b.y);
}
double operator ^(const point &b)const{
return x*b.y-y*b.x;
}
};
struct line{
point s,e;
double k; //斜率
line(){}
line(point _s,point _e):s(_s),e(_e){k=atan2(e.y-s.y,e.x-s.x);}
point operator &(const line &b)const{ //求两直线交点
point res=s;
double t=((s-b.s)^(b.s-b.e))/((s-e)^(b.s-b.e));
res.x+=(e.x-s.x)*t; res.y+=(e.y-s.y)*t;
return res;
}
};
point p[maxn]; //记录最初给的点集
line L[maxn]; //由最初的点集生成直线的集合
point pp[maxn]; //记录半平面交的结果的点集
//半平面交,直线的左边代表有效区域
bool HPIcmp(line a,line b){ //直线排序函数
if(sgn(a.k-b.k)!=) return a.k<b.k; //斜率排序
//斜率相同我也不知道怎么办
return ((a.s-b.s)^(b.e-b.s))<;
}
void HPI(line L[],int n,point res[],int &resn){
//line是半平面交的直线的集合 n是直线的条数 res是结果
//的点集 resn是点集里面点的个数
int tot=n;
sort(L,L+n,HPIcmp);
tot=;
//去掉斜率重复的
for(int i=;i<n;i++) if(sgn(L[i].k-L[i-].k)!=) L[tot++]=L[i];
int head,tail; line Q[maxn]; //双端队列
Q[head=]=L[]; Q[tail=]=L[];
resn=;
for(int i=;i<tot;i++){
if(sgn((Q[tail].e-Q[tail].s)^(Q[tail-].e-Q[tail-].s))==||
sgn((Q[head].e-Q[head].s)^(Q[head+].e-Q[head+].s))==) return;
while(head<tail && (((Q[tail]&Q[tail-])-L[i].s)^(L[i].e-L[i].s))>eps) tail--;
while(head<tail && (((Q[head]&Q[head+])-L[i].s)^(L[i].e-L[i].s))>eps) head++;
Q[++tail]=L[i];
}
while(head<tail && (((Q[tail]&Q[tail-])-Q[head].s)^(Q[head].e-Q[head].s))>eps) tail--;
while(head<tail && (((Q[head]&Q[head-])-Q[tail].s)^(Q[tail].e-Q[tail].e))>eps) head++;
if(tail<=head+) return;
for(int i=head;i<tail; i++) res[resn++]=Q[i]&Q[i+];
if(head<tail-) res[resn++]=Q[head]&Q[tail];
}
double dist(point a,point b){
return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y));
}
void change(point p1,point p2,point &a,point &b,double p){
double len=dist(p1,p2);
double dx=(p1.y-p2.y)*p/len,dy=(p2.x-p1.x)*p/len;
a.x=p1.x+dx;a.y=p1.y+dy;
b.x=p2.x+dx;b.y=p2.y+dy;
}
int main(){
while (cin >> n >> r){
for (int i=;i<n;i++) cin >> p[i].x >> p[i].y;
reverse(p,p+n);
for (int i=;i<n;i++){
point a,b;
change(p[i],p[(i+)%n],a,b,r);
L[i]=line(a,b);
}
int resn;
HPI(L,n,pp,resn); //resn等于0说明移多了
int res1=,res2=;
for (int i=;i<resn;i++)
for (int j=i+;j<resn;j++){
if (dist(pp[i],pp[j])>dist(pp[res1],pp[res2])) res1=i,res2=j;
}
printf("%.5f %.5f %.5f %.5f\n",pp[res1].x,pp[res1].y,pp[res2].x,pp[res2].y);
}
return ;
}

poj3384

5.poj2540 Hotter Colder

传送:http://poj.org/problem?id=2540

题意:在10×10的房间里,有一个点放有宝藏。从(0,0)开始,每次走一个点,告诉你距离宝藏更近or更远,每次输出宝藏可能所在区域的面积。

分析:半平面交解线性规划问题。

 #include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<string>
using namespace std;
const double eps=1e-;
const int maxn=;
int sgn(double x){
if (fabs(x)<eps) return ;
if (x<) return -;
return ;
}
struct point{
double x,y;
point(){}
point(double _x,double _y):x(_x),y(_y){}
point operator +(const point &b)const{
return point(x+b.x,y+b.y);
}
point operator -(const point &b)const{
return point(x-b.x,y-b.y);
}
double operator ^(const point &b)const{
return x*b.y-y*b.x;
}
point operator *(double k)const{
return point(x*k,y*k);
}
};
point pp[maxn];
struct line{
point p,v;
double k;
line(){}
line(point a,point b):p(a),v(b){k=atan2(v.y,v.x);}
bool operator <(const line &b)const{
return k<b.k;
}
};
line L[maxn],q[maxn];
bool left(const line &l,const point &p){
return (l.v^(p-l.p))>; //点在线的上方
}
point pos(const line &a,const line &b){
point x=a.p-b.p;
double t=(b.v^x)/(a.v^b.v);
return a.p+a.v*t;
}
int HalfplaneIntersection(line L[],int n,point poly[]){
sort(L,L+n);
int first,last;
point *p=new point[n];
line *q=new line[n];
q[first=last=] = L[];
for(int i=;i<n;i++){
while(first < last && !left(L[i],p[last-])) last--;
while(first < last && !left(L[i],p[first])) first++;
q[++last] = L[i];
if(fabs((q[last].v^q[last-].v))<eps) {
last--;
if(left(q[last],L[i].p)) q[last] = L[i];
}
if(first<last) p[last-]=pos(q[last-],q[last]);
}
while(first<last && !left(q[first],p[last-])) last--;
if(last-first <=) return ;
p[last]=pos(q[last],q[first]);
int m=;
for(int i=first;i<=last;i++) poly[m++]=p[i];
return m;
}
double dist(point a){
return sqrt(a.x*a.x+a.y*a.y);
}
point normal(point b) {
double len=dist(b);
return point(-b.y/len,b.x/len);
}
double getarea(point p[],int n){
double sum=;
for (int i=;i<n;i++) sum+=(p[i]^(p[(i+)%n]));
return fabs(sum)/;
}
int main(){
L[]=line(point(,),point(,)-point(,));
L[]=line(point(,),point(,)-point(,));
L[]=line(point(,),point(,)-point(,));
L[]=line(point(,),point(,)-point(,));
int cnt=;
double x1=,y1=,x2,y2; string ch;
int f=;
while (cin >> x2 >> y2 >> ch){
point a=point((x1+x2)/,(y1+y2)/),b;
if (!f){cout << "0.00\n";continue;}
if (ch[]=='H') b=normal(point(x1,y1)-point(x2,y2));
else if (ch[]=='C') b=normal(point(x2,y2)-point(x1,y1));
else{cout << "0.00\n";f=;continue;}
L[cnt++]=line(a,b);
int xx=HalfplaneIntersection(L,cnt,pp);
if (xx==){cout << "0.00\n";f=;continue;}
else{
double ans=getarea(pp,xx);
printf("%.2f\n",ans);
if (sgn(ans)==) f=;
}
x1=x2;y1=y2;
}
return ;
}

poj2540

 #include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<string>
using namespace std;
const double eps=1e-;
const int maxn=;
inline double sqr(double x){return x*x;}
int n; double r;
int sgn(double x){
if (fabs(x)<eps) return ;
if (x<) return -;
return ;
}
struct point{
double x,y;
point(){}
point(double _x,double _y):x(_x),y(_y){}
point operator -(const point &b)const{
return point(x-b.x,y-b.y);
}
double operator ^(const point &b)const{
return x*b.y-y*b.x;
}
};
struct line{
point s,e;
double k; //斜率
line(){}
line(point _s,point _e):s(_s),e(_e){k=atan2(e.y-s.y,e.x-s.x);}
point operator &(const line &b)const{ //求两直线交点
point res=s;
double t=((s-b.s)^(b.s-b.e))/((s-e)^(b.s-b.e));
res.x+=(e.x-s.x)*t; res.y+=(e.y-s.y)*t;
return res;
}
};
line L[maxn]; //由最初的点集生成直线的集合
point pp[maxn]; //记录半平面交的结果的点集
//半平面交,直线的左边代表有效区域
bool HPIcmp(line a,line b){ //直线排序函数
if(sgn(a.k-b.k)!=) return a.k<b.k; //斜率排序
return ((a.s-b.s)^(b.e-b.s))<;
}
void HPI(line L[],int n,point res[],int &resn){
int tot=n;
sort(L,L+n,HPIcmp);
tot=;
//去掉斜率重复的
for(int i=;i<n;i++) if(sgn(L[i].k-L[i-].k)!=) L[tot++]=L[i];
int head,tail; line Q[maxn]; //双端队列
Q[head=]=L[]; Q[tail=]=L[];
resn=;
for(int i=;i<tot;i++){
if(sgn((Q[tail].e-Q[tail].s)^(Q[tail-].e-Q[tail-].s))==||
sgn((Q[head].e-Q[head].s)^(Q[head+].e-Q[head+].s))==) return;
while(head<tail && (((Q[tail]&Q[tail-])-L[i].s)^(L[i].e-L[i].s))>eps) tail--;
while(head<tail && (((Q[head]&Q[head+])-L[i].s)^(L[i].e-L[i].s))>eps) head++;
Q[++tail]=L[i];
}
while(head<tail && (((Q[tail]&Q[tail-])-Q[head].s)^(Q[head].e-Q[head].s))>eps) tail--;
while(head<tail && (((Q[head]&Q[head-])-Q[tail].s)^(Q[tail].e-Q[tail].e))>eps) head++;
if(tail<=head+) return;
for(int i=head;i<tail; i++) res[resn++]=Q[i]&Q[i+];
if(head<tail-) res[resn++]=Q[head]&Q[tail];
}
double getarea(point p[],int n){
double sum=;
for (int i=;i<n;i++) sum+=(p[i]^p[(i+)%n]);
return fabs(sum)/;
}
// 将向量 AB 逆时针旋转角度A
inline point Revolve(double cosA, double sinA, point A, point B){
point B_;
B_.x = (B.x - A.x) * cosA - (B.y - A.y) * sinA + A.x;
B_.y = (B.x - A.x) * sinA + (B.y - A.y) * cosA + A.y;
return B_;
}
int main(){
L[]=line(point(,),point(,));
L[]=line(point(,),point(,));
L[]=line(point(,),point(,));
L[]=line(point(,),point(,));
int cnt=;
double x1=,y1=,x2,y2; string ch;
int f=;
while (cin >> x2 >> y2 >> ch){
point a=point((x1+x2)/,(y1+y2)/),b;
int resn;
if (!f){cout << "0.00\n";continue;}
if (ch[]=='C'){
b=Revolve(,,a,point(x2,y2));
}
else if (ch[]=='H'){
b=Revolve(,-,a,point(x2,y2));
}
else{cout << "0.00\n";f=;continue;}
L[cnt++]=line(a,b);
HPI(L,cnt,pp,resn);
if (resn==){cout << "0.00\n";f=;continue;}
else{
double area=getarea(pp,resn);
printf("%.2f\n",area);
}
x1=x2;y1=y2;
}
return ;
}

poj2540(2)