[计算几何]POJ3608 求2个不相交凸包的最短距离

http://acm.pku.edu.cn/JudgeOnline/problem?id=3608

本题要求不相交凸包的最短距离，可以利用旋转卡壳算法来做。首先让2个凸包的点都按逆时针方向排列，然后找出第一个凸包的y坐标最小的点a,第二个凸包的y坐标最大的点b,然后在a,b上做一条与x轴平行的直线，可以证明最短距离只有可能在这2条平行线之间取得。下面的工作就是沿着凸包来旋转这2条平行线了，在旋转的过程中记录最小值。

具体说明在http://cgm.cs.mcgill.ca/~orm/mind2p.html

 

1. Compute the vertex with minimum y coordinate for P (call it yminP) and the vertex with maximum y coordinate for Q (call it ymaxQ).
2. Construct two lines of support LP and LQ for the polygons at yminP and ymaxQ such that the polygons lie to the right of their respective lines of support. Then LP and LQ have opposite direction, and yminP and ymaxQ form an anti-podal pair between the polygons.
3. Compute dist(yminP,ymaxQ) and keep it as the minimum.
4. Rotate the lines clockwise until one of them coincides with an edge of its polygon.
5. If only one line coincides with an edge, then the vertex-edge anti-podal pair distance should be computed along with the new vertex-vertex anti-podal pairA distance. Both distances are compared the current minimum, which is updated if necessary. If both lines of support coincide with edges, then the situation is somewhat more complex. If the edges "overlap", that is if one can construct a line perpendicular to both edges and intersecting both edges (but not at vertices), then the edge-edge distance should be computed. Otherwise the three new vertex-vertex anti-podal pair distances are computed. All distances are compared to the current minimum which is updated if necessary.
6. Repeat steps 4 and 5, until the lines reach (yminP, ymaxQ) again.
7. Output the minimum distance.

在具体实现的过程中，有2点要特别小心：

1.旋转角度的处理。

2.计算2条平行线之间的最短距离最容易出错，要想清楚最短距离到底是什么。

 

//poj 3608
//author: chenkun
//旋转卡壳
#include <stdio.h>
#include <math.h>

const int maxn=10010;
const double eps=1.0e-6;
const double pi=acos(-1.0);
const double inf=1.0e250;

struct Point {
    double x,y;
};
struct Poly {
    int n;
    Point p[maxn];
};
Poly po1,po2;

inline double min(double a,double b) {return a<b?a:b;}
inline double mydist(Point p1,Point p2) {
    return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));
}

double Area(Poly *po) {
    double area=0;
    for(int i=1;i<=po->n;i++) {
        area+=(po->p[i-1].x*po->p[i%po->n].y-po->p[i%po->n].x*po->p[i-1].y);
    }
    return area/2;
}

void changeClockwise(Poly *po) {
    Point tmp;
    for(int i=0;i<=(po->n-1)/2;i++) {
        tmp=po->p[i];
        po->p[i]=po->p[po->n-1-i];
        po->p[po->n-1-i]=tmp;
    }
}

double angle(Point p1,Point p2,double tangle) {
    double ang,tmp;
    Point p;
    p.x=p2.x-p1.x;
    p.y=p2.y-p1.y;
    if(fabs(p.x)<eps) {
        if(p.y>0) ang=pi/2;
        else ang=3*pi/2;
    }else{
        ang=atan(p.y/p.x);
        if(p.x<0) ang+=pi;
    }
    while(ang<0) ang+=2*pi;
    if(ang>=pi) tangle+=pi;
    if(ang>tangle) tmp=ang-tangle;
    else tmp=pi-(tangle-ang);
    while(tmp>=pi) tmp-=pi;
    if(fabs(tmp-pi)<eps) tmp=0;
    return tmp;
}

double disPointToSeg(Point p1,Point p2,Point p3) {
    double a=mydist(p1,p2);
    double b=mydist(p1,p3);
    double c=mydist(p2,p3);
    if(fabs(a+b-c)<eps) return 0;
    if(fabs(a+c-b)<eps||fabs(b+c-a)<eps) return min(a,b);
    double t1=b*b+c*c-a*a;
    double t2=a*a+c*c-b*b;
    if(t1<=0||t2<=0) return min(a,b);

    double area=(p3.x-p1.x)*(p2.y-p1.y)-(p2.x-p1.x)*(p3.y-p1.y);
    return fabs(area)/(c);
}

inline double disPallSeg(Point p1,Point p2,Point p3,Point p4) {
    return min(min(disPointToSeg(p1,p3,p4),disPointToSeg(p2,p3,p4)),min(disPointToSeg(p3,p1,p2),disPointToSeg(p4,p1,p2)));
}

bool init() {
    scanf("%d %d",&po1.n,&po2.n);
    if(po1.n==0&&po2.n==0) return false;
    for(int i=0;i<po1.n;i++) {
        scanf("%lf %lf",&po1.p[i].x,&po1.p[i].y);
    }
    for(int i=0;i<po2.n;i++) {
        scanf("%lf %lf",&po2.p[i].x,&po2.p[i].y);
    }
    return true;
}

int main() {
    //freopen("in.txt","r",stdin);
    while(init()) {
        if(Area(&po1)<0) changeClockwise(&po1);
        if(Area(&po2)<0) changeClockwise(&po2);
        double ymin=inf,ymax=-inf;
        int k1,k2;
        for(int i=0;i<po1.n;i++)
            if(po1.p[i].y<ymin) {
                ymin=po1.p[i].y;
                k1=i;
            }
        for(int i=0;i<po2.n;i++)
            if(po2.p[i].y>ymax) {
                ymax=po2.p[i].y;
                k2=i;
            }
        double langle=0;
        double angle1,angle2;
        double ans=inf;
        double slope=0;
        while(slope<=140) {
            while(langle>=pi) langle-=pi;
            if(fabs(pi-langle)<eps) langle=0;
            angle1=angle(po1.p[k1],po1.p[(k1+1)%po1.n],langle);
            angle2=angle(po2.p[k2],po2.p[(k2+1)%po2.n],langle);
            if(fabs(angle1-angle2)<eps) {
                double d=disPallSeg(po1.p[k1],po1.p[(k1+1)%po1.n],po2.p[k2],po2.p[(k2+1)%po2.n]);
                if(d<ans) ans=d;
                k1=(k1+1)%po1.n;
                k2=(k2+1)%po2.n;
                langle+=angle1;
                slope+=angle1;
            }else if(angle1<angle2) {
                double d=disPointToSeg(po2.p[k2],po1.p[k1],po1.p[(k1+1)%po1.n]);
                if(d<ans) ans=d;
                k1=(k1+1)%po1.n;
                langle+=angle1;
                slope+=angle1;
            }else{
                double d=disPointToSeg(po1.p[k1],po2.p[k2],po2.p[(k2+1)%po2.n]);
                if(d<ans) ans=d;
                k2=(k2+1)%po2.n;
                langle+=angle2;
                slope+=angle2;
            }
        }
        printf("%.5lf\n",ans);
    }
    return 0;
}