【半平面交】 POJ 3384 Feng Shui

先把每条边压缩r的距离，求出半平面交，然后半平面交的最远点对就是答案了。。。要注意最后的点数只有一个时的情况，此时两个圆重合。。但是半平面交求出的平面是不含直线上的点，所以这时半平面交求出的点集为空。。。我的处理方法是压缩每一条边的时候少压缩一点距离，这样子求出的点不会是空集，注意把握好精度就可以了。。。

#include <iostream>  #include <queue>  #include <stack>  #include <map>  #include <set>  #include <bitset>  #include <cstdio>  #include <algorithm>  #include <cstring>  #include <climits>  #include <cstdlib>#include <cmath>#include <time.h>#define maxn 50005#define maxm 3000005#define eps 1e-10#define mod 998244353#define INF 999999999#define lowbit(x) (x&(-x))#define mp mark_pair#define ls o<<1#define rs o<<1 | 1#define lson o<<1, L, mid  #define rson o<<1 | 1, mid+1, R#define debug(x) printf("AA x = %d BB\n", x);//#pragma comment (linker,"/STACK:102400000,102400000")typedef long long LL;//typedef int LL;using namespace std;LL powmod(LL a, LL b){LL res=1,base=a;while(b){if(b%2)res=res*base%mod;base=base*base%mod;b/=2;}return res;}void scanf(int &__x){__x=0;char __ch=getchar();while(__ch==' '||__ch=='\n')__ch=getchar();while(__ch>='0'&&__ch<='9')__x=__x*10+__ch-'0',__ch = getchar();}// headstruct Point{double x, y;Point(double x = 0, double y = 0) : x(x), y(y) { }bool operator < (const Point &a) const {return a.x < x || (a.x == x && a.y < y);}};typedef Point Vector;struct Line{Point P;Vector v;double ang;Line() {}Line(Point P, Vector v) : P(P), v(v) { ang = atan2(v.y, v.x); }bool operator < (const Line &L) const {return ang < L.ang;}};Vector operator + (Vector A, Vector B) {return Vector(A.x + B.x, A.y + B.y);}Vector operator - (Vector A, Vector B){return Vector(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);}int dcmp(double x){if(fabs(x) < eps) return 0;else return x < 0 ? -1 : 1;}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 Cross(Vector A, Vector B){return A.x * B.y - A.y * B.x;}double Area2(Point A, Point B, Point C){return Cross(B - A, C - A);}double PolyonArea(Point *p, int n){double area = 0;for(int i = 1; i < n-1; i++)area += Cross(p[i] - p[0], p[i+1] - p[0]);return area / 2;}Vector Rotate(Vector A, double rad){return Vector(A.x * cos(rad) - A.y * sin(rad), A.x * sin(rad) + A.y * cos(rad));}Vector Normal(Vector A){double L = Length(A);return Vector(-A.y / L, A.x / L);}bool OnLeft(Line L, Point p){return Cross(L.v, p - L.P) > 0;}Point GetIntersection(Line a, Line b){Vector u = a.P - b.P;double t = Cross(b.v, u) / Cross(a.v, b.v);return a.P + a.v * t;}int ConvexHull(Point *p, int n, Point *ch){sort(p, p+n);int m = 0;for(int i = 0; i < n; i++) {while(m > 1 && Cross(ch[m-1] - ch[m-2], p[i] - ch[m-2]) <= 0) m--;ch[m++] = p[i];}int k = m;for(int i = n-2; i >= 0; i--) {while(m > k && Cross(ch[m-1] - ch[m-2], p[i] - ch[m-2]) <= 0) m--;ch[m++] = p[i];}if(n > 1) m--;return m;}Point p[maxn];Line q[maxn];int HalfplaneInersection(Line* L, int n, Point* poly){//sort(L, L+n);int first, last;q[first = last = 0] = L[0];for(int i = 1; i < n; i++) {while(first < last && !OnLeft(L[i], p[last - 1])) last--;while(first < last && !OnLeft(L[i], p[first])) first++;q[++last] = L[i];if(fabs(Cross(q[last].v, q[last-1].v)) < eps) {last--;if(OnLeft(q[last], L[i].P)) q[last] = L[i];}if(first < last) p[last-1] = GetIntersection(q[last-1], q[last]);}while(first < last && !OnLeft(q[first], p[last-1])) last--;if(last - first <= 1) return 0;p[last] = GetIntersection(q[last], q[first]);int m = 0;for(int i = first; i <= last; i++) poly[m++] = p[i];return m;}Point po[maxn], poly[maxn];Vector v[maxn], v2[maxn];Line line[maxn];double r;int n;void read(void){for(int i = n-1; i >= 0; i--) scanf("%lf%lf", &po[i].x, &po[i].y);}void work(void){int ansi, ansj;double mx = 0;for(int i = 0; i < n; i++) {v[i] = po[(i+1)%n] - po[i];v2[i] = Normal(v[i]);}for(int i = 0; i < n; i++) line[i] = Line(po[i] + v2[i] * (r - 0.00001), v[i]);n = HalfplaneInersection(line, n, poly);for(int i = 0; i < n; i++)for(int j = 0; j < n; j++)if(dcmp(Length(poly[i] - poly[j]) - mx) >= 0) {mx = Length(poly[i] - poly[j]);ansi = i, ansj = j;}printf("%.4f %.4f %.4f %.4f\n", poly[ansi].x, poly[ansi].y, poly[ansj].x, poly[ansj].y);}int main(void){while(scanf("%d%lf", &n, &r)!=EOF) {read();work();}return 0;}