题目链接:http://poj.org/problem?id=3384
题意:给定一个凸多边形,用两个半径均为 r 的圆来覆盖凸多边形,求覆盖最大面积时,两圆的圆心坐标,两圆必须完全在多边形内,两圆可以相交。
做法:将凸多边形的每条边内移 r ,求半平面交,得到新的凸多边形,再以距离最远的两个顶点为圆心,这样得到的覆盖面积最大(拿笔画一画就能理解了)。
注意一组数据(discuss中提到的):
4 1
0 0
0 2
2 2
2 0
就是这里错了,RE两次
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
using namespace std;
const double eps = 1e-8;
const int maxn = 105;
int dq[maxn], top, bot, pn, order[maxn], ln;
struct Point {
double x, y;
} p[maxn];
struct Line {
Point a, b;
double angle;
} l[maxn];
int dblcmp(double k) {
if (fabs(k) < eps) return 0;
return k > 0 ? 1 : -1;
}
double multi(Point p0, Point p1, Point p2) {
return (p1.x-p0.x)*(p2.y-p0.y)-(p1.y-p0.y)*(p2.x-p0.x);
}
bool cmp(int u, int v) {
int d = dblcmp(l[u].angle-l[v].angle);
if (!d) return dblcmp(multi(l[u].a, l[v].a, l[v].b)) < 0; //大于0取向量左半部分为半平面,小于0,取右半部分
return d < 0;
}
void getIntersect(Line l1, Line l2, Point& p) {
double dot1,dot2;
dot1 = multi(l2.a, l1.b, l1.a);
dot2 = multi(l1.b, l2.b, l1.a);
p.x = (l2.a.x * dot2 + l2.b.x * dot1) / (dot2 + dot1);
p.y = (l2.a.y * dot2 + l2.b.y * dot1) / (dot2 + dot1);
}
double getDis(Point a, Point b) {
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
bool judge(Line l0, Line l1, Line l2) {
Point p;
getIntersect(l1, l2, p);
return dblcmp(multi(p, l0.a, l0.b)) > 0; //大于小于符号与上面cmp()中注释处相反
}
void addLine(double x1, double y1, double x2, double y2) {
l[ln].a.x = x1; l[ln].a.y = y1;
l[ln].b.x = x2; l[ln].b.y = y2;
l[ln].angle = atan2(y2-y1, x2-x1);
ln++;
}
bool halfPlaneIntersection() {
int i, j;
for (i = 0; i < ln; i++) order[i] = i;
sort(order, order+ln, cmp);
for (i = 1, j = 0; i < ln; i++)
if (dblcmp(l[order[i]].angle-l[order[j]].angle) > 0)
order[++j] = order[i];
ln = j + 1;
dq[0] = order[0];
dq[1] = order[1];
bot = 0;
top = 1;
for (i = 2; i < ln; i++) {
while (bot < top && judge(l[order[i]], l[dq[top-1]], l[dq[top]])) top--;
while (bot < top && judge(l[order[i]], l[dq[bot+1]], l[dq[bot]])) bot++;
dq[++top] = order[i];
}
while (bot < top && judge(l[dq[bot]], l[dq[top-1]], l[dq[top]])) top--;
while (bot < top && judge(l[dq[top]], l[dq[bot+1]], l[dq[bot]])) bot++;
if (bot + 1 >= top) return false;
dq[++top] = dq[bot];
for (pn = 0, i = bot; i < top; i++, pn++)
getIntersect(l[dq[i]], l[dq[i+1]], p[pn]);
return true;
}
/*
将多边形每条边向内推移h
若以向量ab方向,ab右面为半平面,则dx,dy为函数中所写;
若以向量ab方向,ab左面为半平面,则dx,dy为函数中所求dx,dy的相反数
*/
void innerPush(double h) {
double len, dx, dy;
for (int i = 0; i < ln; i++) {
len = getDis(l[i].a, l[i].b);
dx = (l[i].b.y - l[i].a.y) / len * h;
dy = (l[i].a.x - l[i].b.x) / len * h;
l[i].a.x += dx;
l[i].b.x += dx;
l[i].a.y += dy;
l[i].b.y += dy;
}
}
int main()
{
int i, j, n;
double r;
while (scanf ("%d%lf", &n, &r) != EOF) {
for (i = 0; i < n; i++)
scanf ("%lf%lf", &p[i].x, &p[i].y);
for (ln = i = 0; i < n-1; i++)
addLine(p[i].x, p[i].y, p[i+1].x, p[i+1].y);
addLine(p[i].x, p[i].y, p[0].x, p[0].y);
innerPush(r);
halfPlaneIntersection();
double maxl = 0, tmp;
int id1 = 0, id2 = 0;//必须初始化,可能存在为求到的最远两点的距离为0的情况,这里RE了两次
for (i = 0; i < pn; i++) {
for (j = i + 1; j < pn; j++) {
tmp = getDis(p[i], p[j]);
if (dblcmp(tmp-maxl) > 0) {
id1 = i;
id2 = j;
maxl = tmp;
}
}
}
printf ("%.4lf %.4lf %.4lf %.4lf\n", p[id1].x, p[id1].y, p[id2].x, p[id2].y);
}
return 0;
}