POJ 3384 Feng Shui(计算几何の半平面交+最远点对)

Description

Feng shui is the ancient Chinese practice of placement and arrangement of space to achieve harmony with the environment. George has recently got interested in it, and now wants to apply it to his home and bring harmony to it.

There is a practice which says that bare floor is bad for living area since spiritual energy drains through it, so George purchased two similar round-shaped carpets (feng shui says that straight lines and sharp corners must be avoided). Unfortunately, he is unable to cover the floor entirely since the room has shape of a convex polygon. But he still wants to minimize the uncovered area by selecting the best placing for his carpets, and asks you to help.

You need to place two carpets in the room so that the total area covered by both carpets is maximal possible. The carpets may overlap, but they may not be cut or folded (including cutting or folding along the floor border) — feng shui tells to avoid straight lines.

Input

The first line of the input file contains two integer numbers n and r — the number of corners in George’s room (3 ≤ n ≤ 100) and the radius of the carpets (1 ≤ r ≤ 1000, both carpets have the same radius). The following nlines contain two integers xi and yi each — coordinates of the i-th corner (−1000 ≤ xiyi ≤ 1000). Coordinates of all corners are different, and adjacent walls of the room are not collinear. The corners are listed in clockwise order.

Output

Write four numbers x1y1x2y2 to the output file, where (x1y1) and (x2y2) denote the spots where carpet centers should be placed. Coordinates must be precise up to 4 digits after the decimal point.

If there are multiple optimal placements available, return any of them. The input data guarantees that at least one solution exists.

题目大意:给一个凸多边形围成的房子,顺时针给出点,再给两块半径为r的地毯,要求地毯覆盖面积最大且地毯不能切割or折叠,求地毯最大面积覆盖的时候地毯的圆心坐标,任意输出一组解

思路:房子所有边向内移动r,得到一个凸包,凸包上的最远点对即答案之一

PS:数据在http://neerc.ifmo.ru/past/index.html上面有,虽然很多人说这题JPS有问题,但我在WA了无数次之后发现其实还是自己的代码有问题(打错了一个变量名囧)……

暴力枚举最远点对(好牛逼的数据量):

  1 #include <cstdio>
  2 #include <cstring>
  3 #include <cmath>
  4 #include <algorithm>
  5 using namespace std;
  6 
  7 #define EPS 1e-8
  8 #define MAXN 1000
  9 
 10 inline int sgn(double x) {
 11     if(fabs(x) < EPS) return 0;
 12     return x > 0 ? 1 : -1;
 13 }
 14 
 15 struct Point {
 16     double x, y;
 17     Point(double xx = 0, double yy = 0): x(xx), y(yy) {}
 18     bool operator == (const Point &b) const {
 19         return sgn(x - b.x) == 0 && sgn(y - b.y) == 0;
 20     }
 21 };
 22 //cross
 23 inline double operator ^ (const Point &a, const Point &b) {
 24     return a.x * b.y - a.y * b.x;
 25 }
 26 
 27 inline Point operator - (const Point &a, const Point &b) {
 28     return Point(a.x - b.x, a.y - b.y);
 29 }
 30 
 31 struct Line {
 32     Point s, e;
 33     double ag;
 34 };
 35 
 36 struct polygon {
 37     Point v[MAXN];
 38     int n;
 39 } pg, res;
 40 
 41 inline double dist(Point &a, Point &b) {
 42     return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
 43 }
 44 
 45 inline double Cross(Point o, Point s, Point e) {
 46     return (s - o) ^ (e - o);
 47 }
 48 //cross_point
 49 Point operator * (const Line &a, const Line &b) {
 50     Point res;
 51     double u = Cross(a.s, a.e, b.s), v = Cross(a.e, a.s, b.e);
 52     res.x = (b.s.x * v + b.e.x * u)/(u + v);
 53     res.y = (b.s.y * v + b.e.y * u)/(u + v);
 54     return res;
 55 }
 56 
 57 int parallel(Line a, Line b) {
 58     double u = (a.e.x - a.s.x) * (b.e.y - b.s.y) - (a.e.y - a.s.y) * (b.e.x - b.s.x);
 59     return sgn(u) == 0;
 60 }
 61 
 62 inline void set_vector(double x1, double y1, double x2, double y2, Line &v) {
 63     v.s.x = x1; v.s.y = y1;
 64     v.e.x = x2; v.e.y = y2;
 65     v.ag = atan2(y2 - y1, x2 - x1);
 66 }
 67 
 68 Line vct[MAXN], deq[MAXN];
 69 
 70 bool cmp(const Line &a, const Line &b) {
 71     if(sgn(a.ag - b.ag) == 0)
 72         return sgn(Cross(b.s, b.e, a.s)) < 0;
 73     return a.ag < b.ag;
 74 }
 75 
 76 int half_planes_cross(Line *v, int vn) {
 77     int i, n;
 78     //sort(v, v + vn, cmp);
 79     for(i = n = 1; i < vn; ++i) {
 80         if(sgn(v[i].ag - v[i-1].ag) == 0) continue;
 81         v[n++] = v[i];
 82     }
 83     int head = 0, tail = 1;
 84     deq[0] = v[0], deq[1] = v[1];
 85     for(i = 2; i < n; ++i) {
 86         if(parallel(deq[tail - 1], deq[tail]) || parallel(deq[head], deq[head + 1]))
 87             return false;
 88         while(head < tail && sgn(Cross(v[i].s, v[i].e, deq[tail - 1] * deq[tail])) > 0)
 89             --tail;
 90         while(head < tail && sgn(Cross(v[i].s, v[i].e, deq[head] * deq[head + 1])) > 0)
 91             ++head;
 92         deq[++tail] = v[i];
 93     }
 94     while(head < tail && sgn(Cross(deq[head].s, deq[head].e, deq[tail - 1] * deq[tail])) > 0)
 95         --tail;
 96     while(head < tail && sgn(Cross(deq[tail].s, deq[tail].e, deq[head] * deq[head + 1])) > 0)
 97         ++head;
 98     if(tail <= head + 1) return false;
 99     res.n = 0;
100     for(i = head; i < tail; ++i)
101         res.v[res.n++] = deq[i] * deq[i + 1];
102     res.v[res.n++] = deq[head] * deq[tail];
103     res.n = unique(res.v, res.v + res.n) - res.v;
104     res.v[res.n] = res.v[0];
105     return true;
106 }
107 
108 void moving(Line v[], int vn, double r) {
109     for(int i = 0; i < vn; ++i) {
110         double dx = v[i].e.x - v[i].s.x, dy = v[i].e.y - v[i].s.y;
111         dx = dx / dist(v[i].s, v[i].e) * r;
112         dy = dy / dist(v[i].s, v[i].e) * r;
113         v[i].s.x += dy; v[i].e.x += dy;
114         v[i].s.y -= dx; v[i].e.y -= dx;
115     }
116 }
117 
118 int main() {
119     int n;
120     double r;
121     while(scanf("%d%lf", &n, &r) != EOF) {
122         for(int i = 0; i < n; ++i) scanf("%lf%lf", &pg.v[i].x, &pg.v[i].y);
123         pg.v[n] = pg.v[0];
124         for(int i = 0; i < n; ++i)
125             set_vector(pg.v[i].x, pg.v[i].y, pg.v[i+1].x, pg.v[i+1].y, vct[i]);
126         moving(vct, n, r);
127         half_planes_cross(vct, n);
128         int ix = 0, jx = 0;
129         double maxdis = 0;
130         for(int i = 0; i < res.n; ++i) {
131             for(int j = 0; j < res.n; ++j) {
132                 if(i == j) continue;
133                 double t = dist(res.v[i], res.v[j]);
134                 if(sgn(t - maxdis) > 0) {
135                     maxdis = t;
136                     ix = i, jx = j;
137                 }
138             }
139         }
140         printf("%.4f %.4f %.4f %.4f
", res.v[ix].x, res.v[ix].y, res.v[jx].x, res.v[jx].y);
141     }
142 }
View Code

 旋转卡壳求最远点对:

  1 #include <cstdio>
  2 #include <cstring>
  3 #include <cmath>
  4 #include <algorithm>
  5 using namespace std;
  6 
  7 #define EPS 1e-8
  8 #define MAXN 1000
  9 
 10 inline int sgn(double x) {
 11     if(fabs(x) < EPS) return 0;
 12     return x > 0 ? 1 : -1;
 13 }
 14 
 15 struct Point {
 16     double x, y;
 17     Point(double xx = 0, double yy = 0): x(xx), y(yy) {}
 18     bool operator == (const Point &b) const {
 19         return sgn(x - b.x) == 0 && sgn(y - b.y) == 0;
 20     }
 21 };
 22 //cross
 23 inline double operator ^ (const Point &a, const Point &b) {
 24     return a.x * b.y - a.y * b.x;
 25 }
 26 
 27 inline Point operator - (const Point &a, const Point &b) {
 28     return Point(a.x - b.x, a.y - b.y);
 29 }
 30 
 31 struct Line {
 32     Point s, e;
 33     double ag;
 34 };
 35 
 36 struct polygon {
 37     Point v[MAXN];
 38     int n;
 39 } pg, res;
 40 
 41 inline double dist(Point &a, Point &b) {
 42     return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
 43 }
 44 
 45 inline double Cross(Point o, Point s, Point e) {
 46     return (s - o) ^ (e - o);
 47 }
 48 //cross_point
 49 Point operator * (const Line &a, const Line &b) {
 50     Point res;
 51     double u = Cross(a.s, a.e, b.s), v = Cross(a.e, a.s, b.e);
 52     res.x = (b.s.x * v + b.e.x * u)/(u + v);
 53     res.y = (b.s.y * v + b.e.y * u)/(u + v);
 54     return res;
 55 }
 56 
 57 int parallel(Line a, Line b) {
 58     double u = (a.e.x - a.s.x) * (b.e.y - b.s.y) - (a.e.y - a.s.y) * (b.e.x - b.s.x);
 59     return sgn(u) == 0;
 60 }
 61 
 62 inline void set_vector(double x1, double y1, double x2, double y2, Line &v) {
 63     v.s.x = x1; v.s.y = y1;
 64     v.e.x = x2; v.e.y = y2;
 65     v.ag = atan2(y2 - y1, x2 - x1);
 66 }
 67 
 68 Line vct[MAXN], deq[MAXN];
 69 
 70 bool cmp(const Line &a, const Line &b) {
 71     if(sgn(a.ag - b.ag) == 0)
 72         return sgn(Cross(b.s, b.e, a.s)) < 0;
 73     return a.ag < b.ag;
 74 }
 75 
 76 int half_planes_cross(Line *v, int vn) {
 77     int i, n;
 78     //sort(v, v + vn, cmp);
 79     for(i = n = 1; i < vn; ++i) {
 80         if(sgn(v[i].ag - v[i-1].ag) == 0) continue;
 81         v[n++] = v[i];
 82     }
 83     int head = 0, tail = 1;
 84     deq[0] = v[0], deq[1] = v[1];
 85     for(i = 2; i < n; ++i) {
 86         if(parallel(deq[tail - 1], deq[tail]) || parallel(deq[head], deq[head + 1]))
 87             return false;
 88         while(head < tail && sgn(Cross(v[i].s, v[i].e, deq[tail - 1] * deq[tail])) > 0)
 89             --tail;
 90         while(head < tail && sgn(Cross(v[i].s, v[i].e, deq[head] * deq[head + 1])) > 0)
 91             ++head;
 92         deq[++tail] = v[i];
 93     }
 94     while(head < tail && sgn(Cross(deq[head].s, deq[head].e, deq[tail - 1] * deq[tail])) > 0)
 95         --tail;
 96     while(head < tail && sgn(Cross(deq[tail].s, deq[tail].e, deq[head] * deq[head + 1])) > 0)
 97         ++head;
 98     if(tail <= head + 1) return false;
 99     res.n = 0;
100     for(i = head; i < tail; ++i)
101         res.v[res.n++] = deq[i] * deq[i + 1];
102     res.v[res.n++] = deq[head] * deq[tail];
103     res.n = unique(res.v, res.v + res.n) - res.v;
104     res.v[res.n] = res.v[0];
105     return true;
106 }
107 
108 void moving(Line v[], int vn, double r) {
109     for(int i = 0; i < vn; ++i) {
110         double dx = v[i].e.x - v[i].s.x, dy = v[i].e.y - v[i].s.y;
111         dx = dx / dist(v[i].s, v[i].e) * r;
112         dy = dy / dist(v[i].s, v[i].e) * r;
113         v[i].s.x += dy; v[i].e.x += dy;
114         v[i].s.y -= dx; v[i].e.y -= dx;
115     }
116 }
117 
118 int ix, jx;
119 
120 double dia_roataing_calipers() {
121     double dia = 0;
122     ix = jx = 0;
123     int q = 1;
124     for(int i = 0; i < res.n; ++i) {
125         while(sgn(Cross(res.v[i+1], res.v[i], res.v[q+1]) - Cross(res.v[i+1], res.v[i], res.v[q])) > 0)
126             q = (q + 1) % res.n;
127         if(sgn(dist(res.v[i], res.v[q]) - dia) > 0) {
128             dia = dist(res.v[i], res.v[q]);
129             ix = i; jx = q;
130         }
131         if(sgn(dist(res.v[i+1], res.v[q]) - dia) > 0) {
132             dia = dist(res.v[i+1], res.v[q]);
133             ix = i+1; jx = q;
134         }
135     }
136     return dia;
137 }
138 
139 int main() {
140     int n;
141     double r;
142     while(scanf("%d%lf", &n, &r) != EOF) {
143         for(int i = 0; i < n; ++i) scanf("%lf%lf", &pg.v[i].x, &pg.v[i].y);
144         pg.v[n] = pg.v[0];
145         for(int i = 0; i < n; ++i)
146             set_vector(pg.v[i].x, pg.v[i].y, pg.v[i+1].x, pg.v[i+1].y, vct[i]);
147         moving(vct, n, r);
148         half_planes_cross(vct, n);
149         dia_roataing_calipers();
150         printf("%.4f %.4f %.4f %.4f
", res.v[ix].x, res.v[ix].y, res.v[jx].x, res.v[jx].y);
151     }
152 }
View Code
原文地址:https://www.cnblogs.com/oyking/p/3195556.html