旋转卡壳(一)

转载至http://www.cnblogs.com/xdruid/archive/2012/07/01/2572303.html

那么，先提一下最基本最暴力的求凸包直径的方法吧---枚举。。。好吧。。很多问题都可以用 枚举 这个“万能”的方法来解决，过程很简单方便是肯定的，不过在效率上就要差很远了。  要求一个点集的直径，即使先计算出这个点集的凸包，然后再枚举凸包上的点对，这样来求点集直径的话依然会在凸包上点的数量达到O（n）级别是极大的降低它的效率，也浪费了凸包的优美性质。不过在数据量较小或者很适合时，何必要大费周折的用那些麻烦复杂的算法呢，枚举依然是解决问题的很好的方法之一。

 

然后就是今天的旋转卡壳算法了。（那个字念 qia 。。。搞了好久才发现读都读错了。囧。。。）

旋转卡壳可以用于求凸包的直径、宽度，两个不相交凸包间的最大距离和最小距离等。虽然算法的思想不难理解，但是实现起来真的很容易让人“卡壳”。

其实简单来说就是用一对平行线“卡”住凸包进行旋转。

被一对卡壳正好卡住的对应点对称为对踵点，对锺点的具体定义不好说，不过从图上还是比较好理解的。

可以证明对踵点的个数不超过3N/2个 也就是说对踵点的个数是O(N)的

对踵点的个数也是我们下面解决问题时间复杂度的保证。

 

卡壳呢，具体来说有两种情况：

1.

一种是这样，两个平行线正好卡着两个点；

2.

 

一种是这样，分别卡着一条边和一个点。

而第二种情况在实现中比较容易处理，这里就只研究第二种情况。

 

在第二种情况中 我们可以看到 一个对踵点和对应边之间的距离比其他点要大（借用某大神的图··）

 

也就是一个对踵点和对应边所形成的三角形是最大的 下面我们会据此得到对踵点的简化求法。

下面给出一个官方的说明：

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). 

要是真的按这个实现起来就麻烦到吐了。。

根据上面的第二种情况，我们可以得到下面的方法：

如果qa，qb是凸包上最远两点，必然可以分别过qa，qb画出一对平行线。通过旋转这对平行线，我们可以让它和凸包上的一条边重合,如图中蓝色直线，可以注意到，qa是凸包上离p和qb所在直线最远的点。于是我们的思路就是枚举凸包上的所有边，对每一条边找出凸包上离该边最远的顶点，计算这个顶点到该边两个端点的距离，并记录最大的值。

直观上这是一个O(n2)的算法，和直接枚举任意两个顶点一样了。

然而我们可以发现 凸包上的点依次与对应边产生的距离成单峰函数（如下图：）

 

这个性质就很重要啦。

 

根据这个凸包的特性，我们注意到逆时针枚举边的时候，最远点的变化也是逆时针的，这样就可以不用从头计算最远点，而可以紧接着上一次的最远点继续计算。于是我们得到了O(n)的算法。这就是所谓的“旋转”吧！

利用旋转卡壳，我们可以在O（n）的时间内得到凸包的对锺点中的长度最长的点对。

 

又由于最远点对必然属于对踵点对集合 ，那么我们先求出凸包 然后求出对踵点对集合 然后选出距离最大的即可

 

那么具体的代码就很容易实现了，利用叉积，代码只有这么几行的长度：

 

double rotating_calipers(Point *ch,int n) 
{
    int q=1;
    double 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(dist(ch[p],ch[q]),dist(ch[p+1],ch[q+1])));
     }
return ans;
}

其中cross()是计算叉积，因为凸包上距离一条边最远的点和这条边的两个端点构成的三角形面积是最大的。之所以既要更新(ch[p],ch[q])又要更新(ch[p+1],ch[q+1])是为了处理凸包上两条边平行的特殊情况。

下面是POJ 2187, 最基础的模板的题

我的代码:

#include <cstdio>
#include <iostream>
#include <cstring>
#include <cmath>
#include <cstdlib>
#include <algorithm>
#define EPS 1e-8
using namespace std;
const int maxn = 500006;
struct point{
    double x, y;
};
point p[maxn], convex[maxn];
bool cmp(const point p1, const point p2)//比较函数，不用极角，用坐标(x, y); 
{
    return ((p1.y == p2.y && p1.x < p2.x) || p1.y < p2.y);//找最左下角的点为最小点
}
int sgn(double x)
{
    if (fabs(x) < EPS)
        return 0;
    return x < 0 ? -1 : 1;
}
double x_multi(point p1, point p2, point p3)
{
    return ((p3.x - p1.x) * (p2.y - p1.y) - (p2.x - p1.x) * (p3.y - p1.y));
}
double get_distance(point p1, point p2)
{
    return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));
}
double dist2(const point p1, const point p2)//距离的平方
{
    return (p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y);
}
double Max(double a, double b)
{
    return a > b ? a : b;
}
//求凸包, 原来点的集合存在p中, 凸包中的点集合存在convex中, len为凸包中点的个数,n是原来点集的个数
void convex_hull(point *p, point *convex, int n, int &len)
{
    sort(p, p + n, cmp);
    int top = 1;
    convex[0] = p[0];
    convex[1] = p[1];
    //找出凸包的下半部分凸壳 
    for (int i = 2; i < n; i++)
    {
        while (top > 0 && sgn(x_multi(convex[top - 1], convex[top], p[i])) <= 0)//大于0为逆时针,小于0为顺时针
            top--;
        convex[++top] = p[i];
    }
    int tmp = top;
    //找出凸包的上半部分，因为我的比较函数是写的y优先的，所以上下部分，当然也可以以x优先排序，这时候就是左右部分了 
    for (int i = n - 2; i >= 0; i--)
    {
        while (top > tmp && sgn(x_multi(convex[top - 1], convex[top], p[i])) <= 0)//大于0为逆时针,小于0为顺时针
            top--;
        convex[++top] = p[i];//存放凸包中的点
    }
    len = top;
}
//旋转卡壳, 返回最远点对的距离的平方
double rotating_calipers(point *convex, int n)
{
    double ans = 0;
    int q = 1;
    convex[n] = convex[0];
    for (int p = 0; p < n; p++)
    {
        while (x_multi(convex[p], convex[p + 1], convex[q + 1]) > x_multi(convex[p], convex[p + 1], convex[q]))
            q = (q + 1) % n;
        ans = Max(ans, Max(dist2(convex[p], convex[q]), dist2(convex[p + 1], convex[q + 1])));//ans为距离的平方
    }
    return ans;
}
int main()
{
    int n, len;
    while (~scanf("%d", &n))
    {
        for (int i = 0; i < n; i++)
            scanf("%lf %lf", &p[i].x, &p[i].y);
        convex_hull(p, convex, n, len);
        //for (int i = 0; i < len; i++) printf("%.0f, %.0f\n", convex[i].x, convex[i].y);
        printf("%.0f\n", rotating_calipers(convex, len));
    }
    return 0;
}