旋转卡壳（卡尺）算法

申明，本文非笔者原创，原文转载自：http://hi.baidu.com/final_field/item/0c95a87c9299a65d0d0a07f0

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

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
      

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

//计算凸包直径，输入凸包ch，顶点个数为n，按逆时针排列，输出直径的平方
int rotating_calipers(Point *ch,int n)
{
    int q=1,ans=0;
    ch[n]=ch[0];
    for(int p=0;p<n;p++)
     {
        while(cross(ch[p+1],ch[q+1],ch[p])>cross(ch[p+1],ch[q],ch[p]))
            q=(q+1)%n;
        ans=max(ans,max(dist2(ch[p],ch[q]),dist2(ch[p+1],ch[q+1])));            
    }
    return ans; 
}

   很难想象这个看起来那么麻烦的算法只有这么几行代码吧！其中cross函数是计算叉积，可以想成是计算三角形面积，因为凸包上距离一条边最远的点和这条边的两个端点构成的三角形面积是最大的。之所以既要更新(ch[p],ch[q])又要更新(ch[p+1],ch[q+1])是为了处理凸包上两条边平行的特殊情况。

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

poj2187
#include <cmath>
#include <algorithm>
#include <iostream>
using namespace std;
#define MAXN 50005

struct Point
{
    int x, y;
    bool operator < (const Point& _P) const
     {
        return y<_P.y||(y==_P.y&&x<_P.x);
    };
}pset[MAXN],ch[MAXN];

void convex_hull(Point *p,Point *ch,int n,int &len)
{
    sort(p, p+n);
    ch[0]=p[0];
    ch[1]=p[1];
    int top=1;
    for(int i=2;i<n;i++)
     {
        while(top>0&&cross(ch[top],p[i],ch[top-1])<=0)
            top--;
        ch[++top]=p[i];
    }
    int tmp=top;
    for(int i=n-2;i>=0;i--)
     {
        while(top>tmp&&cross(ch[top],p[i],ch[top-1])<=0)
            top--;
        ch[++top]=p[i];
    }
    len=top;    
}

int cross(Point a,Point b,Point o)      
{
    return (a.x - o.x) * (b.y - o.y) - (b.x - o.x) * (a.y - o.y);
}

int dist2(Point a,Point b)
{
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}

int rotating_calipers(Point *ch,int n)
{
    int q=1,ans=0;
    ch[n]=ch[0];
    for(int p=0;p<n;p++)
     {
        while(cross(ch[p+1],ch[q+1],ch[p])>cross(ch[p+1],ch[q],ch[p]))
            q=(q+1)%n;
        ans=max(ans,max(dist2(ch[p],ch[q]),dist2(ch[p+1],ch[q+1])));            
    }
    return ans; 
}

int main()
{
    //freopen("in.txt","r",stdin);
    int n, len;
    while(scanf("%d", &n)!=EOF)
     {
        for(int i = 0;i < n;i++)
         {
            scanf("%d %d",&pset[i].x,&pset[i].y);
        }
        convex_hull(pset,ch,n,len);
        printf("%d\n",rotating_calipers(ch,len));
    }
    return 0;
}

  因相关性关系，有删节。

原文地址：http://www.cppblog.com/staryjy/archive/2010/09/25/101412.html

其他相关的地址：http://blog.youkuaiyun.com/wangyangkobe/archive/2010/12/17/6081975.aspx

                            http://blog.youkuaiyun.com/kaytowin/archive/2010/01/06/5140111.aspx