poj 2187（旋转卡壳）

最新推荐文章于 2020-03-05 20:55:15 发布

qq_30843195

最新推荐文章于 2020-03-05 20:55:15 发布

阅读量241

点赞数

文章标签： poj 2187 旋转卡壳计算几何

旋转卡壳可以用于求凸包的直径、宽度，两个不相交凸包间的最大距离和最小距离等。虽然算法的思想不难理解，但是实现起来真的很容易让人“卡壳”。
   拿凸包直径(也就是凸包上最远的两点的距离）为例，原始的算法是这样子：

Compute the polygon's extreme points in the y direction. Call them ymin and ymax.
Construct two horizontal lines of support through ymin and ymax. Since this is already an anti-podal pair, compute the distance, and keep as maximum.
Rotate the lines until one is flush with an edge of the polygon.
A new anti-podal pair is determined. Compute the new distance, compare to old maximum, and update if necessary.
Repeat steps 3 and 4 until the anti-podal pair considered is (ymin,ymax) again.
Output the pair(s) determining the maximum as the diameter pair(s).

   更具体的可参见http://cgm.cs.mcgill.ca/~orm/rotcal.frame.html


   直接按照这个描述可以实现旋转卡壳算法，但是代码肯定相当冗长。逆向思考，如果q_a，q_b是凸包上最远两点，必然可以分别过q_a，q_b画出一对平行线。通过旋转这对平行线，我们可以让它和凸包上的一条边重合,如图中蓝色直线，可以注意到，q_a是凸包上离p和q_b所在直线最远的点。于是我们的思路就是枚举凸包上的所有边，对每一条边找出凸包上离该边最远的顶点，计算这个顶点到该边两个端点的距离，并记录最大的值。直观上这是一个O(n²)的算法，和直接枚举任意两个顶点一样了。但是注意到当我们逆时针枚举边的时候，最远点的变化也是逆时针的，这样就可以不用从头计算最远点，而可以紧接着上一次的最远点继续计算(详细的证明可以参见上面链接中的论文)。于是我们得到了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])是为了处理凸包上两条边平行的特殊情况。

poj2187要求的是平面点集上的最远点对，实际上就是该点集的凸包的直径。可能该题数据求得的凸包顶点数都不多，所以旋转卡壳算法相比普通的枚举算法并没有明显的优势。完整代码如下。

#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;
    }
}
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()));
    }
     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;
           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());
    }
    //cout << "Hello world!" << endl;
    return 0;
}