您的位置:首页 > 其它

poj3608(旋转卡壳求两凸包间的最短距离)

2017-02-27 22:34 411 查看
/*
translation:
求两个凸包间的最小距离?
solution:
旋转卡壳法。
note:
* 网上给出的资料都差不多,具体如下:
1. 计算凸包P在y轴方向上的最小值记为yminP,和凸包Q在y轴方向上的最大值记为ymaxQ。
2. 建立两条紧贴着yminP, ymaxQ的两条水平的直线LP, LQ。要求他们指向不同的方向。这时候他们就形成了一对anti-podal pair。
3. 计算(yminP,ymaxQ)的距离,并记为minimum.
4. 将两条直线顺时针旋转,直到其中一条遇到凸包的一条边停止。
5. 只要有一条直线遇到了一条边,我们就需要计算新的vertex-vertex anti-podal之间的距离,并同minimum比较,并更新。如果两条直线都分同     一     条     边重合,那么情况就复杂一些了。如果这两条边”重合“,也就是说我们可以画一条垂直于这两条直线的直线并且和凸包上的两个边都相交(不包括在顶点相交),那么我们就需要计算两条直线的距离了。否则的话我们只要计算3个新的顶点到顶点之间的距离。所有的具体都要同minimum比较,并更新minimum 的值。
6. 重复步骤4,5,直到两条直线又回到起始位置为止。
7. 输出最小的距离。

# 之所以在最后旋转卡壳求2次是因为旋转卡壳要求两个直线都已经回到原点后才能结束。但是由于具体代码实现时,结束时不一定保证两条
直线都已经回到原点(因为一个凸包的点数可能较多)。所以可以将两个凸包互换后再处理一次,就能保证正确性。
# 题目所给的点都是按照时针顺序排好的,如果不是的话还需要对其进行排序。
*/
#include <iostream>
#include <cstdio>
#include <vector>
#include <cmath>
#include <algorithm>

using namespace std;
const int maxn = 50000 + 5;
const double INF = 0x3f3f3f3f * 1.0;

struct Point
{
double x, y;
Point(){}
Point(double x_, double y_):x(x_),y(y_){}
} p1[maxn], p2[maxn], t1[maxn], t2[maxn];
typedef Point Vector;
int n, m;

Vector operator + (Vector a, Vector b)  { return Vector(a.x + b.x, a.y + b.y); }
Vector operator - (Point a, Point b)    { return Point(a.x - b.x, a.y - b.y); }
Vector operator * (Vector a, double p)  { return Vector(a.x * p, a.y * p); }
Vector operator / (Vector a, double p)  { return Vector(a.x / p, a.y / p); }

bool operator < (const Point& a, const Point& b)
{
return a.x < b.x || (a.x == b.x && a.y < b.y);
}

const double eps = 1e-10;
int dcmp(double x)
{
if(fabs(x) < eps)    return 0;
else                return x < 0 ? -1 : 1;
}

bool operator == (const Point& a, const Point& b)
{
return dcmp(a.x - b.x) == 0 && dcmp(a.y - b.y) == 0;
}

double dot(Vector a, Vector b)      { return a.x * b.x + a.y * b.y; }
double length(Vector a)             { return sqrt(dot(a, a)); }
double angle(Vector a, Vector b)    { return acos(dot(a, b) / length(a) / length(b)); }
double angle(Vector v)              { return atan2(v.y, v.x); }
double cross(Vector a, Vector b)    { return a.x * b.y - b.x * a.y; }
double dist(Point p1,Point p2)      { return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)); }
double area(Point a, Point b, Point c) { return cross(b - a, c - a); }

double distanceToSegment(Point p, Point a, Point b)
{
if(a.x == b.x && a.y == b.y)	return length(p - a);
Vector v1 = b - a, v2 = p - a, v3 = p - b;

if(dcmp(dot(v1, v2)) < 0)	return length(v2);
else if(dcmp(dot(v1, v3)) > 0)	return length(v3);
else return fabs(cross(v1, v2)) / length(v1);
}

double lineDistance(Point a1, Point b1, Point a2, Point b2)
{
double d1 = min(distanceToSegment(a1, a2, b2), distanceToSegment(b1, a2, b2));
double d2 = min(distanceToSegment(a2, a1, b1), distanceToSegment(b2, a1, b1));
return min(d1, d2);
}

double rotatingCalipers(Point* ch1, Point* ch2, int n, int m)
{
int maxy = 0, miny = 0;
for(int i = 0; i < n; i++)	if(ch1[i].y > ch1[maxy].y)	maxy = i;
for(int i = 0; i < m; i++)	if(ch2[i].y < ch2[miny].y)	miny = i;

ch1
= ch1[0];
ch2[m] = ch2[0];

double ans = INF;
for(int i = 0; i < n; i++) {
double arg;
Vector v = ch1[maxy + 1] - ch1[maxy];
while((arg = cross(v, ch2[miny + 1] - ch1[maxy]) - cross(v, ch2[miny] - ch1[maxy])) > eps)
miny = (miny + 1) % m;
ans = min(ans, lineDistance(ch1[maxy], ch1[maxy + 1], ch2[miny], ch2[miny + 1]));
maxy = (maxy + 1) % n;
}
return ans;
}

int main()
{
//freopen("in.txt", "r", stdin);
while(~scanf("%d%d", &n, &m) && n + m) {
for(int i = 0; i < n; i++)
scanf("%lf%lf", &p1[i].x, &p1[i].y);
for(int j = 0; j < m; j++)
scanf("%lf%lf", &p2[j].x, &p2[j].y);

printf("%.5lf\n", min(rotatingCalipers(p1, p2, n, m), rotatingCalipers(p2, p1, m, n)));
}
return 0;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: