poj 3384 Feng Shui (Half Plane Intersection)-CSDN博客

　　构造半平面交，然后求凸包上最远点对。

　　这题的题意是给出一个凸多边形区域，要求在其中放置两个半径为r的圆（不能超出凸多边形区域），要求求出两个圆心，使得多边形中没有被覆盖的面积最小。反之就是求圆覆盖的区域最大。首先我们可以求出圆心放置的位置的区域，这个要利用半平面交，将原多边形区域向内收缩r的距离。要求两个圆覆盖的区域最大，也就是它们相交的面积最小，也就是两个圆心的距离要尽可能的大。这样就说明了，这题的做法是要求出凸包上面的最远点对。

　　做这题的时候犯了两个错误，一个是没有设置对精度，直接用了cout的默认输出，另一个则是没有想到收缩以后，剩余的多边形的顶点数会少于n。

代码如下：

  1 #include <cstdio>
  2 #include <cstring>
  3 #include <iostream>
  4 #include <algorithm>
  5 #include <vector>
  6 #include <cmath>
  7 
  8 using namespace std;
  9 
 10 struct Point {
 11     double x, y;
 12     Point() {}
 13     Point(double x, double y) : x(x), y(y) {}
 14 } ;
 15 template<class T> T sqr(T x) { return x * x;}
 16 
 17 typedef Point Vec;
 18 Vec operator + (Vec a, Vec b) { return Vec(a.x + b.x, a.y + b.y);}
 19 Vec operator - (Vec a, Vec b) { return Vec(a.x - b.x, a.y - b.y);}
 20 Vec operator * (Vec a, double p) { return Vec(a.x * p, a.y * p);}
 21 Vec operator / (Vec a, double p) { return Vec(a.x / p, a.y / p);}
 22 
 23 const double EPS = 1e-8;
 24 const double PI = acos(-1.0);
 25 inline int sgn(double x) { return (x > EPS) - (x < -EPS);}
 26 bool operator < (Point a, Point b) { return sgn(a.x - b.x) < 0 || sgn(a.x - b.x) == 0 && a.y < b.y;}
 27 bool operator == (Point a, Point b) { return sgn(a.x - b.x) == 0 || sgn(a.y - b.y) == 0;}
 28 
 29 inline double dotDet(Vec a, Vec b) { return a.x * b.x + a.y * b.y;}
 30 inline double crossDet(Vec a, Vec b) { return a.x * b.y - a.y * b.x;}
 31 inline double dotDet(Point o, Point a, Point b) { return dotDet(a - o, b - o);}
 32 inline double crossDet(Point o, Point a, Point b) { return crossDet(a - o, b - o);}
 33 inline double vecLen(Vec x) { return sqrt(dotDet(x, x));}
 34 inline double toRad(double deg) { return deg/ 180.0 * PI;}
 35 inline double angle(Vec v) { return atan2(v.y, v.x);}
 36 Vec normal(Vec x) {
 37     double len = vecLen(x);
 38     return Vec(-x.y, x.x) / len;
 39 }
 40 
 41 struct Poly {
 42     vector<Point> pt;
 43     Poly() { pt.clear();}
 44     ~Poly() {}
 45     Poly(vector<Point> &pt) : pt(pt) {}
 46     Point operator [] (int x) const { return pt[x];}
 47     int size() { return pt.size();}
 48     double area() {
 49         double ret = 0.0;
 50         int sz = pt.size();
 51         for (int i = 1; i < sz; i++) {
 52             ret += crossDet(pt[i], pt[i - 1]);
 53         }
 54         return fabs(ret / 2.0);
 55     }
 56 } ;
 57 
 58 struct DLine {
 59     Point p;
 60     Vec v;
 61     double ang;
 62     DLine() {}
 63     DLine(Point p, Vec v) : p(p), v(v) { ang = atan2(v.y, v.x);}
 64     bool operator < (const DLine &L) const { return ang < L.ang;}
 65     DLine move(double x) {
 66         Vec nor = normal(v);
 67         nor = nor * x;
 68         return DLine(p + nor, v);
 69     }
 70 } ;
 71 
 72 inline bool onLeft(DLine L, Point p) { return crossDet(L.v, p - L.p) > 0;}
 73 Point dLineIntersect(DLine a, DLine b) {
 74     Vec u = a.p - b.p;
 75     double t = crossDet(b.v, u) / crossDet(a.v, b.v);
 76     return a.p + a.v * t;
 77 }
 78 
 79 Poly halfPlane(DLine *L, int n) {
 80     Poly ret = Poly();
 81     sort(L, L + n);
 82     int fi, la;
 83     Point *p = new Point[n];
 84     DLine *q = new DLine[n];
 85     q[fi = la = 0] = L[0];
 86     for (int i = 1; i < n; i++) {
 87         while (fi < la && !onLeft(L[i], p[la - 1])) la--;
 88         while (fi < la && !onLeft(L[i], p[fi])) fi++;
 89         q[++la] = L[i];
 90         if (fabs(crossDet(q[la].v, q[la - 1].v)) < EPS) {
 91             la--;
 92             if (onLeft(q[la], L[i].p)) q[la] = L[i];
 93         }
 94         if (fi < la) p[la - 1] = dLineIntersect(q[la - 1], q[la]);
 95     }
 96     while (fi < la && !onLeft(q[fi], p[la - 1])) la--;
 97     if (la < fi) return ret;
 98     p[la] = dLineIntersect(q[la], q[fi]);
 99     for (int i = fi; i <= la; i++) ret.pt.push_back(p[i]);
100     return ret;
101 }
102 
103 const int N = 111;
104 Point pt[N];
105 DLine dl[N];
106 
107 int main() {
108 //    freopen("in", "r", stdin);
109     int n;
110     double r;
111     while (cin >> n >> r) {
112         for (int i = 0; i < n; i++) {
113             cin >> pt[i].x >> pt[i].y;
114             if (i) dl[i - 1] = DLine(pt[i], pt[i - 1] - pt[i]).move(r + EPS);
115         }
116         dl[n - 1] = DLine(pt[0], pt[n - 1] - pt[0]).move(r + EPS);
117         Poly tmp = halfPlane(dl, n);
118         if (tmp.size() <= 1) {
119             for (int i = 0; i < n; i++) {
120                 if (i) dl[i - 1] = DLine(pt[i], pt[i - 1] - pt[i]).move(r - EPS);
121             }
122             dl[n - 1] = DLine(pt[0], pt[n - 1] - pt[0]).move(r - EPS);
123             tmp = halfPlane(dl, n);
124         }
125         double dis = 0.0;
126         int id[2] = { 0, 0};
127         n = tmp.size();
128         for (int i = 0; i < n; i++) {
129             for (int j = 0; j < n; j++) {
130                 if (dis < vecLen(tmp[i] - tmp[j])) {
131                     dis = vecLen(tmp[i] - tmp[j]);
132                     id[0] = i;
133                     id[1] = j;
134                 }
135             }
136         }
137 //        cout << vecLen(tmp[id[0]] - tmp[id[1]]) << endl;
138         cout.precision(9);
139         cout << tmp[id[0]].x << ' ' << tmp[id[0]].y << ' ' << tmp[id[1]].x << ' ' << tmp[id[1]].y << endl;
140     }
141     return 0;
142 }