POJ - 3608 Bridge Across Islands【旋转卡壳】及一些有趣现象

给两个凸包，求这两个凸包间最短距离

 

旋转卡壳的基础题

因为是初学旋转卡壳，所以找了别人的代码进行观摩。。然而发现很有意思的现象

比如说这个代码（只截取了关键部分）

double solve(Point* P, Point* Q, int n, int m)
{
    int yminP = 0, ymaxQ = 0;
    for (int i = 0; i < n; ++i) if (P[i].y < P[yminP].y) yminP = i;    // P上y坐标最小的顶点
    for (int i = 0; i < m; ++i) if (Q[i].y > Q[ymaxQ].y) ymaxQ = i; // Q上y坐标最大的顶点
    P[n] = P[0];    // 为了方便，避免求余
    Q[m] = Q[0];
    double arg, ans = INF;
    for (int i = 0; i < n; ++i)
    {
        while (arg = cross(P[yminP + 1], Q[ymaxQ + 1], P[yminP]) - cross(P[yminP + 1], Q[ymaxQ], P[yminP]) > EPS) ymaxQ = (ymaxQ + 1) % m;  /*******************/
        if (arg < -EPS) ans = min(ans, point_to_line(P[yminP], P[yminP + 1], Q[ymaxQ]));    // arg!=0不平行
        else ans = min(ans, line_to_line(P[yminP], P[yminP + 1], Q[ymaxQ], Q[ymaxQ + 1]));    // arg=0 平行
        yminP = (yminP + 1) % n;
    }
    return ans;
}
 

我们来看我打星号的那一行代码

while (arg = cross(P[yminP + 1], Q[ymaxQ + 1], P[yminP]) - cross(P[yminP + 1], Q[ymaxQ], P[yminP]) > EPS) ymaxQ = (ymaxQ + 1) % m;

我们知道 “=”，“-”，“>"的优先级是先运算 - 再运算 > 最后运算 =

所以arg的值并不是后边的cross(P[yminP + 1], Q[ymaxQ + 1], P[yminP]) - cross(P[yminP + 1], Q[ymaxQ], P[yminP])

而是后边那个布尔运算的值，也就是说arg只能是0或者1。那么就是说后边arg<-EPS的判断是完全没用的（博主可能想要判断线段Q[ymaxQ] - Q[ymaxQ+1] 和 P[ymaxP] - P[ymaxP+1]是否平行）

我把

if (arg < -EPS) ans = min(ans, point_to_line(P[yminP], P[yminP + 1], Q[ymaxQ]));

和后边的else删掉之后交上去仍然AC。。。

 

我们来看看为什么不会执行那一句话仍然AC，执行这句

ans = min(ans, line_to_line(P[yminP], P[yminP + 1], Q[ymaxQ], Q[ymaxQ + 1]));

显然也包括了上边那句话，只是有更多的冗余判断

所以之前那句话即使不执行仍然可以AC

 

网上很多代码都在犯这个错误，他们是抄的同一份吗ww

 

至于正解，就是求两次旋转卡壳（第二次交换两个凸包的位置），取两个中的最小值

Q: 为什么要执行两次求最小

A: 因为我们进行旋转卡壳时，默认是其中一个凸包的切线正好与某个边重合，然后求另外一凸包离这个边最近的那个点，

因而忽略了这种情况：第一个凸包上的点到第二个凸包某个边距离可能更短。

所以第二次执行时交换两个凸包的位置重新求一次即可。

 

正解程序

#include<cstdio>
#include<cstring>
#include<cmath>
#include<iostream>
#include<algorithm>
#include<set>
#include<map>
#include<stack>
#include<vector>
#include<queue>
#include<string>
#include<sstream>
#define eps 1e-9
#define ALL(x) x.begin(),x.end()
#define INS(x) inserter(x,x.begin())
#define rep(i,j,k) for(int i=j;i<=k;i++)
#define MAXN 20005
#define MAXM 40005
#define INF 0x3fffffff
#define PB push_back
#define MP make_pair
#define X first
#define Y second
#define clr(x,y) memset(x,y,sizeof(x));
using namespace std;
typedef long long LL;
int i,j,k,n,m,x,y,T,ans,big,cas,num,len;
bool flag;

const double pi=acos(-1.0);
int dcmp(double x) {if(fabs(x) < eps) return 0; else return x < 0 ? -1 : 1; }
struct Vector
{
    double x, y;
    Vector (double x=0, double y=0) :x(x),y(y) {}
    Vector operator + (const Vector &B) const { return Vector (x+B.x,y+B.y); }
    Vector operator - (const Vector &B) const { return Vector(x - B.x, y - B.y); }
    Vector operator * (const double &p) const { return Vector(x*p, y*p); }
    Vector operator / (const double &p) const { return Vector(x/p, y/p); }
    double operator * (const Vector &B) const { return x*B.x + y*B.y;}//点积
    double operator ^ (const Vector &B) const { return x*B.y - y*B.x;}//叉积
    bool operator < (const Vector &b) const { return x < b.x || (x == b.x && y < b.y); }
    bool operator ==(const Vector &b) const { return dcmp(x-b.x) == 0 && dcmp(y-b.y) == 0; }
};
typedef Vector Point;
Point Read(){double x, y;scanf("%lf%lf", &x, &y);return Point(x, y);}
double Length(Vector A){ return sqrt(A*A); }//向量的模
double Angle(Vector A, Vector B){return acos(A*B / Length(A) / Length(B)); }//向量的夹角，返回值为弧度
double Area2(Point A, Point B, Point C){ return (B-A)^(C-A); }//向量AB叉乘AC的有向面积
Vector VRotate(Vector A, double rad){return Vector(A.x*cos(rad) - A.y*sin(rad), A.x*sin(rad) + A.y*cos(rad));}//向量A旋转rad弧度
Point  PRotate(Point A, Point B, double rad){return A + VRotate(B-A, rad);}//将B点绕A点旋转rad弧度
Vector Normal(Vector A){double l = Length(A);return Vector(-A.y/l, A.x/l);}//求向量A向左旋转90°的单位法向量,调用前确保A不是零向量

Point GetLineIntersection/*求直线交点，调用前要确保两条直线有唯一交点*/(Point P, Vector v, Point Q, Vector w){double t = (w^(P - Q)) / (v^w);return P + v*t;}//在精度要求极高的情况下，可以自定义分数类
double DistanceToLine/*P点到直线AB的距离*/(Point P, Point A, Point B){Vector v1 = B - A, v2 = P - A;return fabs(v1^v2) / Length(v1);}//不加绝对值是有向距离

double DistanceToSegment(Point P,Point A,Point B){
    if(A == B)  return Length(P-A);
    Vector v1 = B - A,v2 = P - A,v3 = P - B;
    if(dcmp(v1*v2) < 0)       return Length(v2);
    else if(dcmp(v1*v3) > 0)  return Length(v3);
    else    return  fabs((v1^v2))/Length(v1);
}

Point GetLineProjection/*点在直线上的射影*/(Point P, Point A, Point B)
{
    Vector v=B-A;
    return A+v*((v*(P-A))/(v*v));
}

bool OnSegment/*判断点是否在线段上(含端点)*/(Point P,Point a1,Point a2)
{
    Vector v1=a1-P,v2=a2-P;
    if (dcmp(v1^v2)==0 && min(a1.x,a2.x)<=P.x  && P.x<=max(a1.x,a2.x)  && min(a1.y,a2.y)<=P.y && P.y<=max(a1.y,a2.y)) return true;
    return false;
}

bool SegmentInter/*线段相交判定*/(Point a1, Point a2, Point b1, Point b2)
{
    //if (OnSegment(a1,b1,b2) || OnSegment(a2,b1,b2) || OnSegment(b1,a1,a2) || OnSegment(b2,a1,a2)) return 1;
    //如果只判断线段规范相交（不算交点）,上面那句可以删掉
    double c1=(a2-a1)^(b1-a1),c2=(a2-a1)^(b2-a1);
    double c3=(b2-b1)^(a1-b1),c4=(b2-b1)^(a2-b1);
    return dcmp(c1)*dcmp(c2)<0 && dcmp(c3)*dcmp(c4)<0;
}

bool InTri/*判断点是否在三角形内*/(Point P, Point a,Point b,Point c)
{
    if (dcmp(fabs((c-a)^(c-b))-fabs((P-a)^(P-b))-fabs((P-b)^(P-c))-fabs((P-a)^(P-c)))==0) return true;
    return false;
}

double PolygonArea/*求多边形面积，注意凸包P序号从0开始*/(Point *P ,int n)
{
    double ans = 0.0;
    for(int i=1;i<n-1;i++)
        ans+=(P[i]-P[0])^(P[i+1]-P[0]);
    return ans/2;
}
bool CrossOfSegAndLine/*判断线段是否与直线相交*/(Point a1,Point a2,Point b1,Vector b2)
{
    if (OnSegment(b1,a1,a2) || OnSegment(b1+b2,a1,a2)) return true;
    return dcmp(b2^(a1-b1))*dcmp(b2^(a2-b1))<0;
}

Point ch[MAXN];
int ConvexHull(Point* p,int n)
{
    sort(p,p+n);
    int m=0;
    for(int i=0;i<n;++i)
    {
        while(m>1&&((ch[m-1]-ch[m-2])^(p[i]-ch[m-2]))<=0) m--;
        ch[m++]=p[i];
    }
    int k=m;
    for(int i=n-2;i>=0;i--)
    {
        while(m>k&&((ch[m-1]-ch[m-2])^(p[i]-ch[m-2]))<=0) m--;
        ch[m++]=p[i];
    }
    if(n>1) m--;
    for (int i=0;i<m;i++) p[i]=ch[i];
    return m;
}

double Cross(Point A, Point B,Point C)
{
    return (B-A)^(C-A);
}

double dis_pair_seg(Point p1, Point p2, Point p3, Point p4)
{
    return min(min(DistanceToSegment(p1, p3, p4), DistanceToSegment(p2, p3, p4)),
     min(DistanceToSegment(p3, p1, p2), DistanceToSegment(p4, p1, p2)));
}

double rotating_calipers(Point *p1,Point *p2,int n,int m)
{
    int x=0,y=0;
    for (int i=1;i<n;i++) if (p1[i].y<p1[x].y) x=i;
    for (int i=1;i<m;i++) if (p2[i].y>p2[y].y) y=i;

    p1[n]=p1[0];
    p2[m]=p2[0];
    double ans=INF;
    int t=0;
    for (int i=0;i<n;i++)
    {
        while ( (t=dcmp( Cross(p1[x+1],p2[y+1],p1[x]) - Cross(p1[x+1],p2[y],p1[x]))) > 0 )
        {
            y=(y+1)%m;
        }
        if (t<0) ans=min(ans,DistanceToSegment(p2[y],p1[x],p1[x+1]));//不平行
        else ans=min(ans,dis_pair_seg(p1[x],p1[x+1],p2[y],p2[y+1]));//平行
        x=(x+1)%n;
    }
    return ans;
}

Point p1[MAXN],p2[MAXN];
int main()
{
    while (scanf("%d%d",&n,&m),n+m)
    {
        for (i=0;i<n;i++)
        {
            p1[i]=Read();
        }
        n=ConvexHull(p1,n);
        for (i=0;i<m;i++)
        {
            p2[i]=Read();
        }
        m=ConvexHull(p2,m);
        printf("%.9lf\n",min(rotating_calipers(p1,p2,n,m),rotating_calipers(p2,p1,m,n)));  //这里执行两次取最小
    }
}