计算几何（半平面交&&旋转卡壳）

先介绍个写的很好的blog ————Master_Chivu

[Poj 1113] 计算几何之凸包(一) {卷包裹算法}

http://www.cnblogs.com/Booble/archive/2011/02/28/1967179.html

http://www.cnblogs.com/Booble/archive/2011/03/10/1980089.html

[Poj 2187]计算几何之凸包(三) {旋转卡壳初步}

http://www.cnblogs.com/Booble/archive/2011/04/03/2004865.html

旋转卡壳是在上面blog里学到的。

半平面交是在http://apps.hi.baidu.com/share/detail/14889384学到的。

下面是我的模板。

2451（Halfe_Plane_Intersection
ca...POJ的RE提示WA。害我看了两天。以为自己模板又挂了- -。
===================================================================
if ( fabs(seg[i].angle-seg[i-1].angle) > eps )
{

while (deq.size()>=2 && cross(seg[i].sou,seg[i].tar,inter( *deq.begin() , *(deq.begin()+1) )) < -eps)
deq.pop_front();
while (deq.size()>=2 && cross(seg[i].sou,seg[i].tar,inter( *(deq.end()-1),*(deq.end()-2) ) ) < -eps )
deq.pop_back();
deq.push_front(seg[i]);
}
while (deq.size()> 2 && cross((deq.end()-1)->sou,(deq.end()-1)->tar,
inter(*(deq.begin()),*(deq.begin()+1) ) ) < -eps)
deq.pop_front();
while (deq.size()> 2 && cross( deq.begin()->sou, deq.begin()->tar,
inter( *(deq.end()-1), *(deq.end()-2) ) ) < -eps )
deq.pop_back();

一直怀疑划线加粗句是废话，固然- -。去掉后能能毫无压力过掉。
====================================================================
#include<iostream>
#include<cstdio>
#include<deque>
#include<cmath>
#include<algorithm>
#include<cstdlib>
using namespace std;
const int N = 300000;
const int M = 300000;
const double eps = 1e-6;
struct node
{
double x,y;
}poi[N],tmp1,tmp2,zero,tpoi[N];
struct SEG
{
node sou,tar;
double angle;
}seg[M];
deque< SEG > deq;
int Test,n,cnt,flag;
double tmp,s;
double dot(node a,node b)
{
return a.x*b.x+a.y*b.y;
}
double cross(node a,node b,node c)
{
return (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x);
}
node inter( SEG a, SEG b)
{
double s1,s2,ts;
s1 = cross(a.sou,a.tar,b.sou);
s2 = cross(a.sou,a.tar,b.tar);
ts = s1-s2;
if (fabs(ts)<eps) { tmp1.x = 10001; tmp1.y = 10001; return tmp1; }
tmp1.x = b.sou.x * (-s2)/ts + b.tar.x * s1/ts;
tmp1.y = b.sou.y * (-s2)/ts + b.tar.y * s1/ts;
return tmp1;
}
bool cmp(const SEG &a,const SEG &b)
{
if ( fabs(a.angle-b.angle) < eps )
return cross(a.sou,a.tar,b.sou)<0 ;
else
return a.angle<b.angle;

}
void Half_Plane_Intersection()
{
deq.clear();
deq.push_front(seg[0]);
int tt = 1; while (fabs(seg[tt].angle - seg[tt-1].angle)<eps) ++tt;
deq.push_front(seg[tt]);
for (int i = tt+1;i<n;++i)
if ( fabs(seg[i].angle-seg[i-1].angle) > eps )
{

while (deq.size()>=2 && cross(seg[i].sou,seg[i].tar,inter( *deq.begin() , *(deq.begin()+1) )) < -eps)
deq.pop_front();
while (deq.size()>=2 && cross(seg[i].sou,seg[i].tar,inter( *(deq.end()-1),*(deq.end()-2) ) ) < -eps )
deq.pop_back();
deq.push_front(seg[i]);
}
while (deq.size()> 2 && cross((deq.end()-1)->sou,(deq.end()-1)->tar,
inter(*(deq.begin()),*(deq.begin()+1) ) ) < -eps)
deq.pop_front();
while (deq.size()> 2 && cross( deq.begin()->sou, deq.begin()->tar,
inter( *(deq.end()-1), *(deq.end()-2) ) ) < -eps )
deq.pop_back();
cnt = 0;
for ( deque< SEG > ::iterator it = deq.begin(); it!=deq.end(); ++it)
{
deque< SEG > ::iterator itt = it+1;
if (itt == deq.end() ) itt = deq.begin();
tpoi[cnt++] = inter( *it , *itt );
}
s = 0;
for (int i = 0;i<cnt;++i)
s += tpoi[i].x*tpoi[(i+1) % cnt].y - tpoi[(i+1) % cnt].x*tpoi[i].y;
if ( s==0 && cnt<=2 ) puts("0.0");
else printf("%.1lf\n",s/-2 +eps);
}
void work()
{
cin >> n;
for (int i = 0;i<n;++i)
cin >> seg[i].sou.x >> seg[i].sou.y >> seg[i].tar.x >> seg[i].tar.y;
seg[n].sou.x = 0; seg[n].sou.y = 0; seg[n].tar.x = 10000; seg[n].tar.y = 0; ++n;
seg[n].sou.x = 10000; seg[n].sou.y = 0; seg[n].tar.x = 10000; seg[n].tar.y = 10000; ++n;
seg[n].sou.x = 10000; seg[n].sou.y = 10000; seg[n].tar.x = 0; seg[n].tar.y = 10000; ++n;
seg[n].sou.x = 0; seg[n].sou.y = 10000; seg[n].tar.x = 0; seg[n].tar.y = 0 ; ++n;
for (int i =0;i<n;++i)
seg[i].angle = atan2(seg[i].tar.y-seg[i].sou.y,seg[i].tar.x-seg[i].sou.x);
sort(seg,seg+n,cmp);
Half_Plane_Intersection();
}
int main()
{
work();
return 0;
}

下面是旋转卡壳模板

POJ3608 旋转卡壳
给你两个不想交的凸包~求两个凸包间的最短距离。
旋转卡壳模板题。
====================================================================================
int change_angle(int tp1,int tp2,double &ang)
{
double ang1,ang2,a1,a2;
ang1 = atan2( HullA[ (tp1+1) % cnta ].y - HullA[ tp1 ].y , HullA[ (tp1+1) % cnta ].x - HullA[ tp1 ].x);
if (ang1 <0) ang1 = 2*pi + ang1;
a1 = ang1-ang; if (a1<0) a1+=2*pi;
ang2 = atan2( -(HullB[ (tp2+1) % cnta ].y - HullB[ tp2 ].y) , -(HullB[ (tp2+1) % cntb ].x - HullB[ tp2 ].x));
if (ang2 <0) ang2 = 2*pi + ang2;
a2 = ang2-ang; if (a2<0) a2+=2*pi;
if ( fabs( a1 - a2 ) < eps ){ ang = ang1; return 2;}
if ( a1<a2 ) {ang = ang1; return 0;} else {ang = ang2; return 1;}
}
计算平行线下一次跟哪边相贴。
返回1跟A ，返回 2 跟B相贴，返回3，同时相贴。

=====================================================================================
while (flag < (n+m)*5)
{
tmp = change_angle(tp1,tp2,ang);
ans = min( ans , dis(HullB[tp2],HullB[(tp2+1)%cntb],HullA[tp1]) );
ans = min( ans , dis(HullB[tp2],HullB[(tp2+1)%cntb],HullA[(tp1+1) % cnta]) );
ans = min( ans ,dis(HullA[tp1], HullA[(tp1+1) % cnta], HullB[tp2] ) );
ans = min( ans ,dis(HullA[tp1], HullA[(tp1+1) % cnta], HullB[(tp2+1) % cntb] ) );
if (tmp == 0)
tp1 = (tp1+1) % cnta;
else if (tmp == 1)
tp2 = (tp2+1) % cntb;
else
{
tp1 = (tp1+1) % cnta;
tp2 = (tp2+1) % cntb;
}
++flag;
}
旋转卡壳主过程~每次判断当前两条平行线所在凸包边的四个点两两之间距离。 旋转之~
不要偷懒，还是每次都比较四个点比较保险- -。通常不会卡这点常数~
旋转次数嘛~多点就好了。
angle转360就可以了吧。偷懒下。。。
============================================================================================


#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
const double pi = acos(-1.0);
double ans;
const int N = 100000;
const double eps = 1e-8;
int n,cnt,cnta,cntb,tp,m,tmp,tp1,tp2;
struct node
{
double x,y;
} HullA[N],p[N],HullB[N],Hull[N];
double cross(node a,node b, node c)
{
return (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x);
}
bool cmp(const node &a ,const node &b)
{
if (a.y == b.y ) return a.x<b.x;
return a.y<b.y;
}
double dis(node a,node b)
{
return sqrt( (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y) );
}
double dot(node a,node b,node c)
{
return (b.x-a.x)*(c.x-a.x)+(b.y-a.y)*(c.y-a.y);
}
double dis(node a,node b, node c)
{
if (dot(a,b,c)<0 || dot(b,a,c) <0) return min(dis(a,c) , dis(b,c) );
return fabs( cross(a,b,c) ) / dis(a,b);
}
void Calc_Hull(int n)
{
for (int i = 0;i<n;++i)
cin >>p[i].x >> p[i].y;
sort(p,p+n,cmp);
cnt = 0;
for (int i = 0;i<n;++i)
{
while (cnt>1 && cross(Hull[cnt-2],Hull[cnt-1],p[i]) <eps ) cnt--;
Hull[cnt++] = p[i];
}
int t = cnt;
for (int i = n-1;i>=0;i--)
{
while (cnt>t && cross(Hull[cnt-2],Hull[cnt-1],p[i]) <eps) cnt--;
Hull[cnt++] = p[i];
}
--cnt;
}
int change_angle(int tp1,int tp2,double &ang)
{
double ang1,ang2,a1,a2;
ang1 = atan2( HullA[ (tp1+1) % cnta ].y - HullA[ tp1 ].y , HullA[ (tp1+1) % cnta ].x - HullA[ tp1 ].x);
if (ang1 <0) ang1 = 2*pi + ang1;
a1 = ang1-ang; if (a1<0) a1+=2*pi;
ang2 = atan2( -(HullB[ (tp2+1) % cnta ].y - HullB[ tp2 ].y) , -(HullB[ (tp2+1) % cntb ].x - HullB[ tp2 ].x));
if (ang2 <0) ang2 = 2*pi + ang2;
a2 = ang2-ang; if (a2<0) a2+=2*pi;
if ( fabs( a1 - a2 ) < eps ){ ang = ang1; return 2;}
if ( a1<a2 ) {ang = ang1; return 0;} else {ang = ang2; return 1;}
}
int main()
{
int Test = 0;
while (cin >> n >> m)
{
if (n == 0) break;
Calc_Hull(n);
cnta = cnt;
for (int i = 0;i<cnta;++i) HullA[i] = Hull[i];
Calc_Hull(m);
cntb = cnt;
for (int i = 0;i<cntb;++i) HullB[i] = Hull[i];
tp1 = 0;tp2 = 0;
for (int i = 1;i<cntb;++i) if ( !cmp(HullB[i],HullB[tp2]) ) tp2 = i;
ans = dis( HullA[0], HullB[tp2] );
int flag = 0;
double ang = 0;
while (flag < (n+m)*5)
{
tmp = change_angle(tp1,tp2,ang);
ans = min( ans , dis(HullB[tp2],HullB[(tp2+1)%cntb],HullA[tp1]) );
ans = min( ans , dis(HullB[tp2],HullB[(tp2+1)%cntb],HullA[(tp1+1) % cnta]) );
ans = min( ans ,dis(HullA[tp1], HullA[(tp1+1) % cnta], HullB[tp2] ) );
ans = min( ans ,dis(HullA[tp1], HullA[(tp1+1) % cnta], HullB[(tp2+1) % cntb] ) );
if (tmp == 0)
tp1 = (tp1+1) % cnta;
else if (tmp == 1)
tp2 = (tp2+1) % cntb;
else
{
tp1 = (tp1+1) % cnta;
tp2 = (tp2+1) % cntb;
}
++flag;
}
printf("%.5lf\n",ans);
}
return 0;
}
posted @ 2011-04-09 10:34  Whimsy  阅读(530)  评论(0)  编辑  收藏  举报