更具体的可参见http://cgm.cs.mcgill.ca/~orm/rotcal.frame.html
直接按照这个描述可以实现旋转卡壳算法,但是代码肯定相当冗长。逆向思考,如果qa,qb是凸包上最远两点,必然可以分别过qa,qb画出一对平行线。通过旋转这对平行线,我们可以让它和凸包上的一条边重合,如图中蓝色直线,可以注意到,qa是凸包上离p和qb所在直线最远的点。于是我们的思路就是枚举凸包上的所有边,对每一条边找出凸包上离该边最远的顶点,计算这个顶点到该边两个端点的距离,并记录最大的值。直观上这是一个O(n2)的算法,和直接枚举任意两个顶点一样了。但是注意到当我们逆时针枚举边的时候,最远点的变化也是逆时针的,这样就可以不用从头计算最远点,而可以紧接着上一次的最远点继续计算(详细的证明可以参见上面链接中的论文)。于是我们得到了O(n)的算法。
int rotating_calipers() //卡壳
{
int q=1;
int ans=0;
stack[top]=0;
for(int i=0;i<top;i++)
{
while( (p[stack[i+1]]-p[stack[i]])*(p[stack[q+1]]-p[stack[i]])>(p[stack[i+1]]-p[stack[i]])*(p[stack[q]]-p[stack[i]]))
q=(q+1)%(top);
ans=max(ans , max((p[stack[i]]-p[stack[q]]).len2(),(p[stack[i+1]]-p[stack[q+1]]).len2()));
}
return ans;
}
很难想象这个看起来那么麻烦的算法只有这么几行代码吧!其中cross函数是计算叉积,可以想成是计算三角形面积,因为凸包上距离一条边最远的点和这条边的两个端点构成的三角形面积是最大的。之所以既要更新(ch[p],ch[q])又要更新(ch[p+1],ch[q+1])是为了处理凸包上两条边平行的特殊情况。
poj 2187 要求的是平面点集上的最远点对,实际上就是该点集的凸包的直径。可能该题数据求得的凸包顶点数都不多,所以旋转卡壳算法相比普通的枚举算法并没有明显的优势。完整代码如下。
#include<iostream>
#include<cmath>
#define PI acos(-1.0)
const double eps = 1e-6;
int sgn(double a) {
return (a>eps) - (a<-eps);
}
struct point{
int x,y;
point (int xx=0,int yy=0){ //构造函数
x=xx; y=yy;
}
point operator -(point b){
return point(x-b.x,y-b.y);
}
int operator *(point b){
return (x*b.y-y*b.x);
}
int len2(){
return x*x+y*y;
}
void input(){
scanf("%d%d",&x,&y);
}
}p[50005];
int stack[50005],top,n;
int cmp(const void *a,const void *b) //逆时针排序 返回正数要交换
{
struct point *c=(struct point *)a;
struct point *d=(struct point *)b;
int k=(*c-p[0])*(*d-p[0]);
if(k<0) return 1;
else if(k==0 && ((*c-p[0]).len2()-(*d-p[0]).len2())>=0)
return 1;
else return -1;
}
void graham(int n) //形成凸包
{
int i;
for(i=0;i<=1;i++)
stack[i]=i;
top=1;
for(i=2;i<n;i++)
{
while(top>0 && (p[stack[top]]-p[stack[top-1]])*(p[i]-p[stack[top-1]])<=0)
top--;
stack[++top]=i;
}
int tmp=top;
for(i=n-2;i>=0;i--)
{
while(top>tmp && (p[stack[top]]-p[stack[top-1]])*(p[i]-p[stack[top-1]])<=0)
top--;
stack[++top]=i;
}
/* for(i=0;i<n;i++)
printf("**%d %d\n",p[i].x,p[i].y);
printf("\n%d\n",top);
for(i=0;i<top;i++)
printf("**%d %d\n",p[stack[i]].x,p[stack[i]].y);*/
}
int max(int a,int b)
{
return a>b?a:b;
}
int rotating_calipers() //卡壳
{
int q=1;
int ans=0;
stack[top]=0;
for(int i=0;i<top;i++)
{
while( (p[stack[i+1]]-p[stack[i]])*(p[stack[q+1]]-p[stack[i]])>
(p[stack[i+1]]-p[stack[i]])*(p[stack[q]]-p[stack[i]]))
q=(q+1)%(top);
ans=max(ans , max((p[stack[i]]-p[stack[q]]).len2(),(p[stack[i+1]]-p[stack[q+1]]).len2()));
// printf("%d stack[i]=%d stack[q]=%d %d %d %d\n",i,stack[i],stack[q],
// (p[stack[i]]-p[stack[q]]).len2(),(p[stack[i+1]]-p[stack[q+1]]).len2(),ans);
}
return ans;
}
int main()
{
int i;
while(scanf("%d",&n)!=EOF)
{
for(i=0;i<n;i++)
p[i].input();
int u=0;
for(i=1;i<n;i++) //找左下角的点
if(p[i].y<p[u].y || (p[i].y==p[u].y&&p[i].x<p[u].x) )
u=i;
//swap
point tmp=p[0];
p[0]=p[u];
p[u]=tmp;
qsort(p+1,n-1,sizeof(p[0]),cmp);
graham(n);
printf("%d\n",rotating_calipers());
}
return 0;
}