计算几何：Most Distant Point from the Sea(半平面交)

001:Most Distant Point from the Sea
查看提交统计提问
总时间限制: 5000ms 内存限制: 65536kB
描述
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.

输入
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 (xi, yi)–(xi+1, yi+1) (1 ≤ i ≤ n − 1) and the line segment (xn, yn)–(x1, y1) 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.

输出
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−5). You may output any number of digits after the decimal point, provided that the above accuracy condition is satisfied.

样例输入
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
样例输出
5000.000000
494.233641
34.542948
0.353553

法向量平移，二分查找验证，半平面交

#include <cmath>
#include <iomanip>
#include <iostream>
#include <string>
#include <utility>
#include <vector>
using namespace std;
#define inff 0x3f3f3f3f
int n;
double INF = 1e20;
double EPS = 1e-6; //精度不是越高越好
bool IsZero(double x) { return -EPS < x && x < EPS; }
int Sign(double x)
{ // 判断 x 是大于0,等于0还是小于0
    return fabs(x) < EPS ? 0 : x > 0 ? 1 : -1;
}
struct CPoint
{
    double x;
    double y;
    CPoint() {}
    CPoint(double xx, double yy)
    {
        x = xx;
        y = yy;
    }
    void init(double xx, double yy)
    {
        x = xx;
        y = yy;
    }
    bool operator==(CPoint &a)
    {
        return IsZero(x - a.x) && IsZero(y - a.y);
    }
};
struct CVector
{
    double x;
    double y;
    CVector() {}
    CVector(double xx, double yy)
    {
        x = xx;
        y = yy;
    }
    CVector operator*(double m)
    {
        x = m * x;
        y = m * y;
        return *this;
    }
    CVector(CPoint &a, CPoint &b)
    {
        x = a.x - b.x;
        y = a.y - b.y;
    }
};
struct CLine
{
    CPoint a;
    CPoint b;
    CLine() {}
    void init(CPoint xx, CPoint yy)
    {
        a = xx;
        b = yy;
    }
    CLine(CPoint xx, CPoint yy)
    {
        a = xx;
        b = yy;
    }
};
CVector operator-(CPoint b, CPoint a) { return CVector(b.x - a.x, b.y - a.y); } // c = a – b;
CVector operator*(double k, CVector p) { return CVector(k * p.x, k * p.y); }
double operator^(CVector p, CVector q) { return p.x * q.y - q.x * p.y; }
CVector operator/(CVector p, double q) { return CVector(p.x / q, p.y / q); }
double area(CVector p, CVector q) { return (p ^ q) / 2; }
double PI = acos(-1);
double operator*(CVector p, CVector q) { return p.x * q.x + p.y * q.y; }
double length(CVector p) { return sqrt(p * p); }
double dist(CPoint p, CPoint q) { return length(p - q); }
CPoint operator+(CPoint a, CVector p) { return CPoint(a.x + p.x, a.y + p.y); }
double dist(CPoint p, CLine l) { return fabs((p - l.a) ^ (l.b - l.a)) / length(l.b - l.a); }
CLine p, q;
CPoint intersect(CLine l, CLine m)
{
    double x = area(m.a - l.a, l.b - l.a);
    double y = area(l.b - l.a, m.b - l.a);
    if (IsZero(x + y))
    {
        CPoint ss(1001, -1001);
        if (IsZero(dist(l.a, m)))
        {
            ss.y -= 1;
            return ss;
        }
        else
        {
            return ss;
        }
    }
    return m.a + ((x / (x + y)) * (m.b - m.a));
}
double Cross(const CVector &v1, const CVector &v2)
{
    //叉积
    return v1.x * v2.y - v1.y * v2.x;
}
bool IsFormalCross(CPoint p1, CPoint p2, CPoint p3, CPoint p4)
{ //判断p1->p2和p3->p4是否规范相交
    return Cross(p2 - p1, p3 - p1) * Cross(p2 - p1, p4 - p1) < -EPS && Cross(p4 - p3, p1 - p3) * Cross(p4 - p3, p2 - p3) < -EPS;
}
double upx[25], upy[25], downx[25], downy[25];
CPoint ups[25], downs[25];

struct Seg //线段
{
    CPoint a, b; //向量是 a->b ,即 b-a
    Seg(const CPoint &aa, const CPoint &bb) : a(aa), b(bb)
    {
    }
    //直线两点式方程 (y-y1)/(y2-y1) = (x-x1)/(x2-x1)
    double getX(double y)
    { //给定y坐标，求直线上的 x坐标
        return (y - a.y) / (b.y - a.y) * (b.x - a.x) + a.x;
    }
    double getY(double x)
    { //给定x坐标，求直线上的y坐标
        return (x - a.x) / (b.x - a.x) * (b.y - a.y) + a.y;
    }
};
bool FLessEq(double a, double b)
{
    return b - a > EPS || IsZero(b - a);
}
bool PointInSeg(CPoint p, Seg L)
{
    double tmp = Cross(L.a - p, L.a - L.b);
    if (!IsZero(tmp))
        return false;
    if (FLessEq(min(L.a.x, L.b.x), p.x) && FLessEq(p.x, max(L.a.x, L.b.x)) && FLessEq(min(L.a.y, L.b.y), p.y) && FLessEq(p.y, max(L.a.y, L.b.y)))
        return true;
    return false;
}
double Distance(const CPoint &p1, const CPoint &p2)
{
    return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));
}

pair<int, CPoint> CrossPoint(Seg s1, Seg s2)
{
    CPoint p1 = s1.a;
    CPoint p2 = s1.b;
    CPoint p3 = s2.a;
    CPoint p4 = s2.b;
    double a1 = area(p3 - p1, p4 - p1);
    double a2 = area(p4 - p2, p3 - p2); //这些顺序不能乱
    if (Cross(p2 - p1, p3 - p1) * Cross(p2 - p1, p4 - p1) < -EPS && Cross(p4 - p3, p1 - p3) * Cross(p4 - p3, p2 - p3) < -EPS)
    { // 规范相交
        return make_pair(0, p1 + (p2 - p1) * (a1 / (a1 + a2)));
    }
    if (!(IsZero(Cross(p2 - p1, p3 - p4))))
    {                             //不平行,不共线
        if (p1 == p3 || p1 == p4) //端点重合且不平行，不共线
            return make_pair(1, p1);
        if (p2 == p3 || p2 == p4)
            return make_pair(1, p2);
        if (PointInSeg(p1, s2))
            return make_pair(2, p1);
        if (PointInSeg(p2, s2))
            return make_pair(2, p2);
        if (PointInSeg(p3, s1))
            return make_pair(3, p3);
        if (PointInSeg(p4, s1))
            return make_pair(3, p4);
        CPoint crossPoint = p1 + (p2 - p1) * (a1 / (a1 + a2));
        if (PointInSeg(crossPoint, s1))
            return make_pair(8, crossPoint);
        if (PointInSeg(crossPoint, s2))
            return make_pair(9, crossPoint);
        return make_pair(4, crossPoint);
        // 直线和线段也无交点，不平行，不共线，两直线交点是second
    }
    CLine ss(s2.a, s2.b);
    if (!IsZero(dist(p1, ss)))
        return make_pair(5, CPoint(0, 0)); //平行
    //下面都是共线，且有公共点 if (PointInSeg(p1, s2)) return make_pair(6, p1);
    if (PointInSeg(p2, s2))
        return make_pair(6, p2);
    if (PointInSeg(p3, s1))
        return make_pair(6, p3);
    if (PointInSeg(p4, s1))
        return make_pair(6, p4);
    return make_pair(7, CPoint(0, 0)); //共线，且无公共点
}
typedef vector<CPoint> Polygon;
int CutPolygon(const Polygon &src, CPoint a, CPoint b, Polygon &result)
{ //用直线a->b 切割src,返回其左侧的多边形，放到result //src里的点什么顺序无所谓
    int n = src.size();
    result.clear();
    for (int i = 0; i < n; ++i)
    {
        CPoint c = src[i];
        CPoint d = src[(i + 1) % n];
        if (Sign(Cross(b - a, c - a)) >= 0) //叉积>=0
            result.push_back(c);
        pair<int, CPoint> r = CrossPoint(Seg(c, d), Seg(a, b));
        //c,d在 a->b上，或者 c->d重合于 a->b这两种情况都不要考 虑，因为会在上面两行处理了
        if (r.first == 0 || r.first == 8 || r.first == 3)
            //规范相交或虽非规范相交但是交点在c->d上
            result.push_back(r.second);
    }
    return result.size();
}
Polygon src;
Polygon result[105];
bool is(double mid)
{
    int ss = src.size();
    result[0] = src;
    for (int i = 0; i < ss; i++)
    {
        int j = i;
        CVector h(src[(j + 1)%ss], src[j]);
        double mmmm = h.x;
        h.x = -h.y;
        h.y = mmmm;
        CPoint x = src[j] + h / length(h) * mid;
        CPoint y = src[(j + 1)%ss] + h / length(h) * mid;
        if(CutPolygon(result[j], x, y, result[(j+1)%ss]) == 0)
        {
            return false;
        }
    }
    return true;
}



int main()
{
  
    CPoint tmp;
    double tx, ty;
    double right, left, mid;
    while (cin >> n && n)
    {
        src.clear();
        for (int i = 0; i < n; i++)
        {
            cin >> tx >> ty;
            tmp.init(tx, ty);
            src.push_back(tmp);
            result[i].clear();
        }
        left = 0;
        right = 10000;
        while(right - 0.000001 > left)
        {
           // cout << mid << endl;
            mid = (left + right) / 2;
            if(is(mid))
            {
                left = mid;
            }
            else 
            {
                right = mid - 0.000001;
            }
        }
        printf("%.6lf\n", left);
    }
}