【PKU3525】半平面交二分求多边形核

Most Distant Point from the Sea

Time Limit: 5000MS		Memory Limit: 65536K
Total Submissions: 3985		Accepted: 1862		Special Judge

Description

The main land of Japan called Honshu is an island surrounded by the sea. In such an island, it is natural to ask a question: “Where is the most distant point from the sea?” The answer to this question for Honshu was found in 1996. The most distant point is located in former Usuda Town, Nagano Prefecture, whose distance from the sea is 114.86 km.

In this problem, you are asked to write a program which, given a map of an island, finds the most distant point from the sea in the island, and reports its distance from the sea. In order to simplify the problem, we only consider maps representable by convex polygons.

Input

The input consists of multiple datasets. Each dataset represents a map of an island, which is a convex polygon. The format of a dataset is as follows.

n
x1		y1
	⋮
xn		yn

Every input item in a dataset is a non-negative integer. Two input items in a line are separated by a space.

n in the first line is the number of vertices of the polygon, satisfying 3 ≤ n ≤ 100. Subsequent n lines are the x- and y-coordinates of the n vertices. Line segments (x_i, y_i)–(x_i+1, y_i+1) (1 ≤ i ≤ n − 1) and the line segment (x_n, y_n)–(x₁, y₁) form the border of the polygon in counterclockwise order. That is, these line segments see the inside of the polygon in the left of their directions. All coordinate values are between 0 and 10000, inclusive.

You can assume that the polygon is simple, that is, its border never crosses or touches itself. As stated above, the given polygon is always a convex one.

The last dataset is followed by a line containing a single zero.

Output

For each dataset in the input, one line containing the distance of the most distant point from the sea should be output. An output line should not contain extra characters such as spaces. The answer should not have an error greater than 0.00001 (10⁻⁵). You may output any number of digits after the decimal point, provided that the above accuracy condition is satisfied.

Sample Input

4
0 0
10000 0
10000 10000
0 10000
3
0 0
10000 0
7000 1000
6
0 40
100 20
250 40
250 70
100 90
0 70
3
0 0
10000 10000
5000 5001
0

Sample Output

5000.000000
494.233641
34.542948
0.353553

Source

Japan 2007

求一个多边形内的最大圆半径~

二分枚举多边形各个边向内移动的距离r，如果存在核，即存在cut出来的多边形，那么说明存在一个半径为r的圆可以放在多边形里

#define DeBUG
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <cstdlib>
#include <algorithm>
#include <vector>
#include <stack>
#include <queue>
#include <string>
#include <set>
#include <sstream>
#include <map>
#include <list>
#include <bitset>
using namespace std ;
#define zero {0}
#define INF 0x3f3f3f3f
#define EPS 1e-10
#define TRUE true
#define FALSE false
typedef long long LL;
const double PI = acos(-1.0);
//#pragma comment(linker, "/STACK:102400000,102400000")
inline int sgn(double x)
{
    return fabs(x) < EPS ? 0 : (x < 0 ? -1 : 1);
}
#define N 100005
// Point class
struct Point;
typedef Point Vec;
struct Point
{
    double x, y;
    Point () {}
    Point(double a, double b)
    {
        x = a;
        y = b;
    }
};
Vec operator + (const Vec &a, const Vec &b) //点加法
{
    return Vec(a.x + b.x, a.y + b.y);
}
Vec operator - (const Vec &a, const Vec &b) //点减法
{
    return Vec(a.x - b.x, a.y - b.y);
}
Vec operator * (const Vec &a, const double &p) //点与常数相乘
{
    return Vec(a.x * p, a.y * p);
}
bool operator == (const Vec &a, const Point &b) //点相等判断
{
    return sgn(a.x - b.x) == 0 && sgn(a.y - b.y) == 0;
}
bool operator < (const Vec &a, const Point &b) //平面直角坐标系中左下方的为小
{
    return a.x < b.x || (sgn(a.x - b.x) == 0 && a.y < b.y);
}
inline double crossDet(Vec a, Vec b)//叉乘
{
    return a.x * b.y - a.y * b.x;
}
inline double crossDet(Point o, Point a, Point b)//向量叉乘
{
    return crossDet(a - o, b - o);
}
inline double dotDet(Vec a, Vec b)//点乘
{
    return a.x * b.x + a.y * b.y;
}
inline double vecLen(Vec a)//求一点到原点距离
{
    return sqrt(dotDet(a, a));
}
Vec normal(Vec a)//逆时针转90度的一个a的单位向量
{
    double len = vecLen(a);
    return Vec(- a.y / len, a.x / len);
}
// Polygon class
struct Poly//平面类
{
    vector<Point> pt;//保存平面对应的端点
    Poly()
    {
        pt.clear();
    }
    ~Poly() {}
    Poly(vector<Point> pt0): pt(pt0) {} //使用vector进行初始化
    Point operator [](int x) const//重载[]返回对应端点
    {
        return pt[x];
    }
    int size()//返回对应平面点数
    {
        return pt.size();
    }
    double area()//得到面积,凹凸多边形均可
    {
        double ret = 0.0;
        int sz = pt.size();
        for (int i = 0; i < sz; i++)
        {
            int j = (i + 1) % sz;
            ret += crossDet(pt[i], pt[j]);
        }
        return fabs(ret / 2.0);
    }
};
struct DLine
{
    Point p;
    Vec v;
    double ang;
    DLine() {}
    DLine(Point p, Vec v): p(p), v(v)
    {
        ang = atan2(v.y, v.x);
    }
    bool operator < (const DLine &L) const
    {
        return ang < L.ang;
    }
    DLine move(double x)   // while x > 0 move to v's left
    {
        Vec nor = normal(v);
        nor = nor * x;
        return DLine(p + nor, v);
    }
} ;
inline bool onLeft(DLine L, Point p)
{
    return crossDet(L.v, p - L.p) > 0;
}
Point dLineIntersect(DLine a, DLine b)
{
    Vec u = a.p - b.p;
    double t = crossDet(b.v, u) / crossDet(a.v, b.v);
    return a.p + a.v * t;
}

Poly halfPlane(DLine *L, int n)
{
    Poly ret = Poly();
    sort(L, L + n);
    int fi, la;
    Point *p = new Point[n];
    DLine *q = new DLine[n];
    q[fi = la = 0] = L[0];
    for (int i = 1; i < n; i++)
    {
        while (fi < la && !onLeft(L[i], p[la - 1])) la--;
        while (fi < la && !onLeft(L[i], p[fi])) fi++;
        q[++la] = L[i];
        if (fabs(crossDet(q[la].v, q[la - 1].v)) < EPS)
        {
            la--;
            if (onLeft(q[la], L[i].p)) q[la] = L[i];
        }
        if (fi < la) p[la - 1] = dLineIntersect(q[la - 1], q[la]);
    }
    while (fi < la && !onLeft(q[fi], p[la - 1])) la--;
    if (la - fi <= 1) return ret;
    p[la] = dLineIntersect(q[la], q[fi]);
    for (int i = fi; i <= la; i++) ret.pt.push_back(p[i]);
    return ret;
}
DLine L[205];
Point p[205];
int n;
bool solve(double r)
{
    DLine ll[205];
    // printf("%.3lf\n", r);
    for (int i = 0; i < n; i++)
    {
        ll[i] = L[i].move(r);//向有向线段左边移动r个单位
    }
    Poly poly = halfPlane(ll, n);
    if (poly.size() > 0)return 1;
    else
        return 0;
}
int main()
{
#ifdef DeBUGs
    freopen("C:\\Users\\Sky\\Desktop\\1.in", "r", stdin);
#endif

    while (scanf("%d", &n),n)
    {
        for (int i = 0; i < n; i++)
        {
            scanf("%lf%lf", &p[i].x, &p[i].y);
        }
        p[n] = p[0];
        for (int i = 0; i < n; i++)
        {
            L[i] = DLine(p[i], p[i + 1] - p[i]);
        }
        double l = 0.0;
        double r = 100000.0;
        double mid;
        while (l + EPS <= r)
        {
            mid = (l + r) / 2.0;
            if (solve(mid))
                l = mid;
            else
                r = mid;
        }
        printf("%.10f\n", l);
    }

    return 0;
}