POJ 3384 Feng Shui 半平面交

#include <iostream>  
#include <cstdio>  
#include <cstring>  
#include <string>  
#include <algorithm>  
#include <cstdlib>  
#include <cmath>  
#include <map>  
#include <sstream>  
#include <queue>  
#include <vector>  
#define MAXN 111111  
#define MAXM 211111  
#define PI acos(-1.0)  
#define eps 1e-8  
#define INF 1000000001  
using namespace std;  
int dblcmp(double d)  
{  
    if (fabs(d) < eps) return 0;  
    return d > eps ? 1 : -1;  
}  
struct point  
{  
    double x, y;  
    point(){}  
    point(double _x, double _y):  
    x(_x), y(_y){};  
    void input()  
    {  
        scanf("%lf%lf",&x, &y);  
    }  
    double dot(point p)  
    {  
        return x * p.x + y * p.y;  
    }  
    double distance(point p)  
    {  
        return hypot(x - p.x, y - p.y);  
    }  
    point sub(point p)  
    {  
        return point(x - p.x, y - p.y);  
    }  
    double det(point p)  
    {  
        return x * p.y - y * p.x;  
    }  
    bool operator == (point a)const  
    {  
        return dblcmp(a.x - x) == 0 && dblcmp(a.y - y) == 0;  
    }  
    bool operator < (point a)const  
    {  
        return dblcmp(a.x - x) == 0 ? dblcmp(y - a.y) < 0 : x < a.x;  
    }  
  
}p[MAXN];  
struct line  
{  
    point a,b;  
    line(){}  
    line(point _a,point _b)  
    {  
        a=_a;  
        b=_b;  
    }  
    bool parallel(line v)  
    {  
        return dblcmp(b.sub(a).det(v.b.sub(v.a))) == 0;  
    }  
    point crosspoint(line v)  
    {  
        double a1 = v.b.sub(v.a).det(a.sub(v.a));  
        double a2 = v.b.sub(v.a).det(b.sub(v.a));  
        return point((a.x * a2 - b.x * a1) / (a2 - a1), (a.y * a2 - b.y * a1) / (a2 - a1));  
    }  
    bool operator == (line v)const  
    {  
        return (a == v.a) && (b == v.b);  
    }  
};  
struct halfplane:public line  
{  
    double angle;  
    halfplane(){}  
    //表示向量 a->b逆时针(左侧)的半平面  
    halfplane(point _a, point _b)  
    {  
        a = _a;  
        b = _b;  
    }  
    halfplane(line v)  
    {  
        a = v.a;  
        b = v.b;  
    }  
    void calcangle()  
    {  
        angle = atan2(b.y - a.y, b.x - a.x);  
    }  
    bool operator <(const halfplane &b)const  
    {  
        return angle < b.angle;  
    }  
};  
struct polygon  
{  
    int n;  
    point p[MAXN];  
    line l[MAXN];  
    double area;  
    void getline()  
    {  
        for (int i = 0; i < n; i++)  
        {  
            l[i] = line(p[i], p[(i + 1) % n]);  
        }  
    }  
    void getarea()  
    {  
        area = 0;  
        int a = 1, b = 2;  
        while(b <= n - 1)  
        {  
            area += p[a].sub(p[0]).det(p[b].sub(p[0]));  
            a++;  
            b++;  
        }  
        area = fabs(area) / 2;  
    }  
}convex;  
struct halfplanes  
{  
    int n;  
    halfplane hp[MAXN];  
    point p[MAXN];  
    int que[MAXN];  
    int st, ed;  
    void push(halfplane tmp)  
    {  
        hp[n++] = tmp;  
    }  
    void unique()  
    {  
        int m = 1, i;  
        for (i = 1; i < n;i++)  
        {  
            if (dblcmp(hp[i].angle - hp[i - 1].angle))hp[m++] = hp[i];  
            else if (dblcmp(hp[m - 1].b.sub(hp[m - 1].a).det(hp[i].a.sub(hp[m - 1].a)) > 0))hp[m - 1] = hp[i];  
        }  
        n = m;  
    }  
    bool halfplaneinsert()  
    {  
        int i;  
        for (i = 0; i < n; i++) hp[i].calcangle();  
        sort(hp, hp + n);  
        unique();  
        que[st = 0] = 0;  
        que[ed = 1] = 1;  
        p[1] = hp[0].crosspoint(hp[1]);  
        for (i = 2; i < n; i++)  
        {  
            while (st < ed && dblcmp((hp[i].b.sub(hp[i].a).det(p[ed].sub(hp[i].a)))) < 0) ed--;  
            while (st < ed && dblcmp((hp[i].b.sub(hp[i].a).det(p[st + 1].sub(hp[i].a)))) < 0) st++;  
            que[++ed] = i;  
            if (hp[i].parallel(hp[que[ed - 1]])) return false;  
            p[ed] = hp[i].crosspoint(hp[que[ed - 1]]);  
        }  
        while (st < ed && dblcmp(hp[que[st]].b.sub(hp[que[st]].a).det(p[ed].sub(hp[que[st]].a))) < 0) ed--;  
        while (st < ed && dblcmp(hp[que[ed]].b.sub(hp[que[ed]].a).det(p[st + 1].sub(hp[que[ed]].a))) < 0) st++;  
        if (st + 1 >= ed)return false;  
        return true;  
    }  
    void getconvex(polygon &con)  
    {  
        p[st] = hp[que[st]].crosspoint(hp[que[ed]]);  
        con.n = ed - st + 1;  
        int j = st, i = 0;  
        for (; j <= ed; i++, j++)  
        {  
            con.p[i] = p[j];  
        }  
    }  
}h;  
int T;  
int n;  
line getmove(point a, point b, double mid)  
{  
    double x = a.x - b.x;  
    double y = a.y - b.y;  
    double L = a.distance(b);  
    point ta = point(mid * y / L + a.x, a.y - mid * x / L);  
    point tb = point(mid * y / L + b.x, b.y - mid * x / L);  
    return line(ta, tb);  
}  
double r;  
int main()  
{  
    int cas = 0;  
    while(scanf("%d%lf", &n, &r) != EOF)  
    {  
        for(int i = 0; i < n; i++) p[i].input();  
        h.n = 0;  
  
        for(int i = 0; i < n; i++)  
        {  
            line tmp = getmove(p[(i + 1) % n], p[i], r);  
            h.push(halfplane(tmp));  
        }  
        h.halfplaneinsert();  
        h.getconvex(convex);  
        int id1 = 0, id2 = 0;  
        double mx = 0;  
        for(int i = 0; i < convex.n; i++)  
            for(int j = i + 1; j < convex.n; j++)  
            {  
                double len = convex.p[i].distance(convex.p[j]);  
                if(dblcmp(len - mx) > 0) id1 = i, id2 = j, mx = len;  
            }  
        printf("%f %f %f %f\n", convex.p[id1].x, convex.p[id1].y, convex.p[id2].x, convex.p[id2].y);  
    }  
    return 0;  
}