http://poj.org/problem?id=3608
给定两个非连接(比如不相交)的凸多边形 P 和 Q, 目标是找到拥有最小距离的点对 (p,q) (p 属于 P 且 q 属于 Q)。
事实上, 多边形非连接十分重要, 因为我们所说的多边形包含其内部。 如果多边形相交, 那么最小距离就变得没有意义了。 然而, 这个问题的另一个版本, 凸多边形顶点对间最小距离对于相交和非相交的情况都有解存在。
回到我们的主问题: 直观的, 确定最小距离的点不可能包含在多边形的内部。 与最大距离问题相似, 我们有如下结论:
两个凸多边形 P 和 Q 之间的最小距离由多边形间的对踵点对确立。 存在凸多边形间的三种多边形间的对踵点对, 因此就有三种可能存在的最小距离模式:
“顶点-顶点”的情况
“顶点-边”的情况
“边-边”的情况
换句话说, 确定最小距离的点对不一定必须是顶点。 下面的三个图例表明了以上结论:
给定结果, 一个基于旋转卡壳的算法自然而然的产生了:
考虑如下的算法, 算法的输入是两个分别有 m 和 n 个顺时针给定顶点的凸多边形 P 和 Q。
计算 P 上 y 坐标值最小的顶点(称为 yminP ) 和 Q 上 y 坐标值最大的顶点(称为 ymaxQ)。
为多边形在 yminP 和 ymaxQ 处构造两条切线 LP 和 LQ 使得他们对应的多边形位于他们的右侧。 此时 LP 和 LQ 拥有不同的方向, 并且 yminP 和 ymaxQ 成为了多边形间的一个对踵点对。
计算距离(yminP,ymaxQ) 并且将其维护为当前最小值。
顺时针同时旋转平行线直到其中一个与其所在的多边形的边重合。
如果只有一条线与边重合, 那么只需要计算“顶点-边”对踵点对和“顶点-顶点”对踵点对距离。 都将他们与当前最小值比较, 如果小于当前最小值则进行替换更新。 如果两条切线都与边重合, 那么情况就更加复杂了。 如果边“交叠”, 也就是可以构造一条与两条边都相交的公垂线(但不是在顶点处相交), 那么就计算“边-边”距离。 否则计算三个新的“顶点-顶点”对踵点对距离。 所有的这些距离都与当前最小值进行比较, 若小于当前最小值则更新替换。
重复执行步骤4和步骤5, 直到新的点对为(yminP,ymaxQ)。
输出最大距离。
旋转卡壳模式保证了所有的对踵点对(和所有可能的子情况)都被考虑到。 此外, 整个算法拥有现行的时间复杂度, 因为(除了初始化), 只有与顶点数同数量级的操作步数需要执行。
最小距离和最大距离的问题表明了旋转卡壳模型可以用在不同的条件下(与先前的直径和宽度问题比较)。 这个模型可以应用于两个多边形的问题中。
“最小盒子”问题(最小面积外接矩形)通过同一多边形上两个正交切线集合展示了另一种条件下旋转卡壳的应用。
#include <cmath>
#include <cstdio>
#include <string.h>
#include <algorithm>
#include <iostream>
typedef double TYPE;
const double Epsilon=1.0e-6;
const double pi=acos(-1.0);
const double inf=1.0e250;
//空间中的点,可以用来作为二维点来用
struct POINT {/*验证*/
TYPE x; TYPE y; TYPE z;
//需要什么变量就将什么变量添加进去
void read() { scanf ("%lf %lf", &x, &y); }
//同上
void write(){ printf("%f %f\n", x, y); }
POINT() : x(0), y(0), z(0) {};
POINT(TYPE _x_, TYPE _y_, TYPE _z_ = 0)
: x(_x_), y(_y_), z(_z_) {};
//要用 G++ 提交 ,可以不用这个
POINT operator =(const POINT &A){
x = A.x;
y = A.y;
z = A.z;
}
POINT operator-(POINT &a)
{
POINT t;
t.x = x - a.x;
t.y = y - a.y;
t.z = z - a.z;
return t;
}
POINT operator+(POINT &a)
{
POINT t;
t.x = x + a.x;
t.y = y + a.y;
t.z = z + a.y;
return t;
}
};
// 多边形 ,逆时针或顺时针给出x,y
struct POLY {/*验证*/
//n个点
int n;
//x,y为点的指针,首尾必须重合
POINT *p;
POLY() : n(0), p(NULL){};
POLY(int _n_){
n = _n_;
p = new POINT[n + 1];
}
POLY(int _n_, const POINT * _p_ ) {
n = _n_;
p = new POINT[n + 1];
memcpy(p, _p_, n*sizeof(POINT));
p[n] = _p_[0];
}
};
inline double Min(double a,double b) {return a<b?a:b;}
// planar points' distance
// 两个点的距离
TYPE Distance(const POINT & a, const POINT & b) { /*验证*/
return
sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y)
* (a.y - b.y) + (a.z - b.z) * (a.z - b.z));
}
//求多边形面积 拍好序的点
TYPE Area(const POLY & poly) { /*验证*/
if ( poly.n < 3)
return TYPE(0);
double s = poly.p[0].y * (poly.p[poly.n - 1].x - poly.p[1].x);
for (int i = 1; i < poly.n; i++) {
s += poly.p[i].y * (poly.p[i - 1].x - poly.p[(i + 1) % poly.n].x);
}
return s/2;
}
/*
//顺时针为 负
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)<Epsilon) {
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)<Epsilon) tmp=0;
return tmp;
}
double disPointToSeg(POINT p1,POINT p2,POINT p3) {
double a=Distance(p1,p2);
double b=Distance(p1,p3);
double c=Distance(p2,p3);
if(fabs(a+b-c)<Epsilon) return 0;
if(fabs(a+c-b)<Epsilon||fabs(b+c-a)<Epsilon) 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)));
}
void hhhh(POLY po1, POLY po2) {
//面积为负的话这说明点排列为顺时针
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)<Epsilon) 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)<Epsilon) {
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("%.5f\n",ans);
}
int main(){
int n1,n2;
while( scanf("%d %d",&n1,&n2)) {
POLY po1(n1);
POLY po2(n2);
if(n1==0&&n2==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);
}
hhhh(po1, po2);
}
}