强烈推荐此题。半平面交算法的一个应用。
具体做法是,把多边形的每条边向内平移r单位长度,用这些线段所在直线和原多边形作半平面交,得到的区域就是半径为r的圆放入多边形的可行域。可以证明这个区域一定是凸的,或者退化为一条线段,或一个点。那么,我们就可以在这个区域上求最远点对啦。
我的做法是O(n
2)的。应该存在O(nlogn)的做法,因为都是凸多边形,每次半平面交只有最多两个交点,可二分,而最后的求最远点对可以旋转卡壳。比赛的时候时间少,就写了个暴力O(n
2)的。
/**//*************************************************************************
Author: WHU_GCC
Created Time: 2007-9-23 12:02:01
File Name: pku3384.cpp
Description:
************************************************************************/
#include <iostream>
#include <cmath>
using namespace std;
#define out(x) (cout<<#x<<": "<<x<<endl)
const int maxint=0x7FFFFFFF;
typedef long long int64;
const int64 maxint64 = 0x7FFFFFFFFFFFFFFFLL;
template<class T>void show(T a, int n){for(int i=0; i<n; ++i) cout<<a[i]<<' '; cout<<endl;}
template<class T>void show(T a, int r, int l){for(int i=0; i<r; ++i)show(a[i],l);cout<<endl;}
const double eps = 1e-10;
const int maxn = 200;
struct point
{
double x, y;
};
struct cp
{
int n;
point p[maxn];
};
double dist(point a, point b)
{
return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
point intersectL(double a1, double b1, double c1, double a2, double b2, double c2)
{
point ret;
ret.y = (a1 * c2 - c1 * a2) / (b1 * a2 - a1 * b2);
if (fabs(a2) < eps)
ret.x = -(b1 * ret.y + c1) / a1;
else
ret.x = -(b2 * ret.y + c2) / a2;
return ret;
}
bool isEqual(point inpA, point inpB)
{
return (fabs(inpA.x - inpB.x) < eps && fabs(inpA.y - inpB.y) < eps);
}
double Cross(point inpA, point inpB, point inpC)
{
return (inpB.x - inpA.x) * (inpC.y - inpA.y) - (inpC.x - inpA.x) * (inpB.y - inpA.y);
}
void get_line(point inpA, point inpB, double &a1, double &b1, double &c1)
{
a1 = inpB.y - inpA.y;
b1 = inpA.x - inpB.x;
c1 = inpA.y * (inpB.x - inpA.x) - inpA.x * (inpB.y - inpA.y);
}
cp cut(point inpA, point inpB, cp incp)
{
cp ret;
point cross;
int i, j;
double t1, t2;
double a1, b1, c1, a2, b2, c2;
ret.n = 0;
for (i = 0; i < incp.n; i++)
{
j = i + 1;
t1 = Cross(inpA, inpB, incp.p[i]);
t2 = Cross(inpA, inpB, incp.p[j]);
if (t1 < eps && t2 < eps)
{
ret.p[ret.n++] = incp.p[i];
ret.p[ret.n++] = incp.p[j];
}
else if (t1 > eps && t2 > eps)
continue;
else
{
get_line(inpA, inpB, a1, b1, c1);
get_line(incp.p[i], incp.p[j], a2, b2, c2);
cross = intersectL(a1, b1, c1, a2, b2, c2);
if (t1 < eps)
{
ret.p[ret.n++] = incp.p[i];
ret.p[ret.n++] = cross;
}
else
{
ret.p[ret.n++] = cross;
ret.p[ret.n++] = incp.p[j];
}
}
}
if (ret.n == 0)
return ret;
for (i = 1, j = 1; i < ret.n; i++)
if (!isEqual(ret.p[i - 1], ret.p[i]))
ret.p[j++] = ret.p[i];
ret.n = j;
if (ret.n != 1 && isEqual(ret.p[ret.n - 1], ret.p[0]))
ret.n--;
ret.p[ret.n] = ret.p[0];
return ret;
}
int main()
{
int n, r;
cp input, ret;
scanf("%d%d", &n, &r);
input.n = n;
for (int i = 0; i < n; i++)
scanf("%lf%lf", &input.p[i].x, &input.p[i].y);
input.p[n] = input.p[0];
ret = input;
for (int i = 0; i < n; i++)
{
point ta, tb, tt;
tt.x = input.p[i + 1].y - input.p[i].y;
tt.y = input.p[i].x - input.p[i + 1].x;
double k = r / sqrt(tt.x * tt.x + tt.y * tt.y);
tt.x = tt.x * k;
tt.y = tt.y * k;
ta.x = input.p[i].x + tt.x;
ta.y = input.p[i].y + tt.y;
tb.x = input.p[i + 1].x + tt.x;
tb.y = input.p[i + 1].y + tt.y;
ret = cut(ta, tb, ret);
}
double ans = -1;
double ans_x1, ans_y1, ans_x2, ans_y2;
for (int i = 0; i < ret.n; i++)
for (int j = 0; j < ret.n; j++)
{
double t = dist(ret.p[i], ret.p[j]);
if (t > ans)
{
ans = t;
ans_x1 = ret.p[i].x;
ans_y1 = ret.p[i].y;
ans_x2 = ret.p[j].x;
ans_y2 = ret.p[j].y;
}
}
printf("%.4lf %.4lf %.4lf %.4lf\n", ans_x1, ans_y1, ans_x2, ans_y2);
return 0;
}