旋转卡壳

旋转卡壳

关注
发信
关注(28)粉丝(399)

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

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

旋转卡壳的基础题

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

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

double solve(Point* P, Point* Q, int n, int m)

{

    int yminP = , ymaxQ = ;

    for (int i = ; i < n; ++i) if (P[i].y < P[yminP].y) yminP = i;    // P上y坐标最小的顶点

    for (int i = ; i < m; ++i) if (Q[i].y > Q[ymaxQ].y) ymaxQ = i; // Q上y坐标最大的顶点

    P[n] = P[];    // 为了方便，避免求余

    Q[m] = Q[];

    double arg, ans = INF;

    for (int i = ; i < n; ++i)

    {

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

        if (arg < -EPS) ans = min(ans, point_to_line(P[yminP], P[yminP + ], Q[ymaxQ]));    // arg!=0不平行

        else ans = min(ans, line_to_line(P[yminP], P[yminP + ], Q[ymaxQ], Q[ymaxQ + ]));    // arg=0 平行

        yminP = (yminP + ) % n;

    }

    return ans;

}

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

while (arg = cross(P[yminP + ], Q[ymaxQ + ], P[yminP]) - cross(P[yminP + ], Q[ymaxQ], P[yminP]) > EPS) ymaxQ = (ymaxQ + ) % 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 + ], Q[ymaxQ]));

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

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

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

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

所以之前那句话即使不执行仍然可以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 ; else return x <  ? - : ; }

struct Vector

{

    double x, y;

    Vector (double x=, double y=) :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) ==  && dcmp(y-b.y) == ; }

};

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) < )       return Length(v2);

    else if(dcmp(v1*v3) > )  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)== && 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)< && dcmp(c3)*dcmp(c4)<;

}

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)))==) return true;

    return false;

}

double PolygonArea/*求多边形面积，注意凸包P序号从0开始*/(Point *P ,int n)

{

    double ans = 0.0;

    for(int i=;i<n-;i++)

        ans+=(P[i]-P[])^(P[i+]-P[]);

    return ans/;

}

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

}

Point ch[MAXN];

int ConvexHull(Point* p,int n)

{

    sort(p,p+n);

    int m=;

    for(int i=;i<n;++i)

    {

        while(m>&&((ch[m-]-ch[m-])^(p[i]-ch[m-]))<=) m--;

        ch[m++]=p[i];

    }

    int k=m;

    for(int i=n-;i>=;i--)

    {

        while(m>k&&((ch[m-]-ch[m-])^(p[i]-ch[m-]))<=) m--;

        ch[m++]=p[i];

    }

    if(n>) m--;

    for (int i=;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=,y=;

    for (int i=;i<n;i++) if (p1[i].y<p1[x].y) x=i;

    for (int i=;i<m;i++) if (p2[i].y>p2[y].y) y=i;

    p1[n]=p1[];

    p2[m]=p2[];

    double ans=INF;

    int t=;

    for (int i=;i<n;i++)

    {

        while ( (t=dcmp( Cross(p1[x+],p2[y+],p1[x]) - Cross(p1[x+],p2[y],p1[x]))) >  )

        {

            y=(y+)%m;

        }

        if (t<) ans=min(ans,DistanceToSegment(p2[y],p1[x],p1[x+]));//不平行

        else ans=min(ans,dis_pair_seg(p1[x],p1[x+],p2[y],p2[y+]));//平行

        x=(x+)%n;

    }

    return ans;

}

Point p1[MAXN],p2[MAXN];

int main()

{

    while (scanf("%d%d",&n,&m),n+m)

    {

        for (i=;i<n;i++)

        {

            p1[i]=Read();

        }

        n=ConvexHull(p1,n);

        for (i=;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)));  //这里执行两次取最小

    }

}
05-06 21:27
Powered by LMLPHP ©2024 旋转卡壳 0.078707