UVALive 3890 Most Distant Point from the Sea(凸包最大内接圆)

最新推荐文章于 2023-09-02 15:53:36 发布

原创最新推荐文章于 2023-09-02 15:53:36 发布 · 1.1k 阅读

0 ·

CC 4.0 BY-SA版权

计算几何专栏收录该内容

19 篇文章

订阅专栏

本文介绍了一种求解凸多边形中最大内切圆半径的方法，通过半平面交技术和二分查找算法实现。该方法首先确定凸包上的每条边对应的半平面，并沿垂直单位法向量平移预设距离，然后检查平移后的半平面交集是否为空来逐步逼近最大内切圆的半径。

一个n个点的凸多边形，求多边形中离多边形边界最远的距离。实际上就是求凸包最大内接圆的半径。

利用半平面交求解，每次二分枚举半径d，然后将凸包每条边所代表的半平面沿其垂直单位法向量平移d，看所有平移后的半平面的交集是否为空。

#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<fstream>
#include<sstream>
#include<bitset>
#include<vector>
#include<string>
#include<cstdio>
#include<cmath>
#include<stack>
#include<queue>
#include<stack>
#include<map>
#include<set>
#define FF(i, a, b) for(int i=a; i<b; i++)
#define FD(i, a, b) for(int i=a; i>=b; i--)
#define REP(i, n) for(int i=0; i<n; i++)
#define CLR(a, b) memset(a, b, sizeof(a))
#define debug puts("**debug**")
#define LL long long
#define PB push_back
#define eps 1e-10
using namespace std;

struct Point
{
    double x, y;
    Point (double x=0, double y=0):x(x), y(y) {}
};
typedef Point Vector;

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); }

bool operator < (const Point& a, const Point& b)
{
    return a.x < b.x || (a.x == b.x && a.y < b.y);
}

int dcmp(double x)
{
    if(fabs(x) < eps) return 0;
    return x < 0 ? -1 : 1;
}

bool operator == (const Point& a, const Point& b)
{
    return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0;
}

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 Angel(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(Vector A, Vector B, Vector C) { return Cross(B-A, C-A); }
Vector Normal(Vector a) //a向量的垂直法向量
{
     return Vector(-a.y/Length(a), a.x/Length(a));
}

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;
    }
};

//点p在半平面的左边
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;
}

const int maxn = 200;
Point p[maxn], poly[maxn];
Line L[maxn];
Vector v[maxn], v2[maxn];
int n;

//半平面交
Point pp[maxn];
Line qq[maxn];
int HalfplaneIntersection(Line* L, int n, Point* poly)
{
    sort(L, L+n);
    int first, last;

    qq[first=last=0] = L[0];
    FF(i, 1, n)
    {
        while(first < last && !onLeft(L[i], pp[last-1])) last--;
        while(first < last && !onLeft(L[i], pp[first])) first++;
        qq[++last] = L[i];
        if(fabs(Cross(qq[last].v, qq[last-1].v)) < eps)
        {
            last--;
            if(onLeft(qq[last], L[i].p)) qq[last] = L[i];
        }
        if(first < last) pp[last-1] = GetIntersection(qq[last-1], qq[last]);
    }
    while(first < last && !onLeft(qq[first], pp[last-1])) last--;
    if(last-first <= 1) return 0;
    pp[last] = GetIntersection(qq[last], qq[first]);

    int m = 0;
    FF(i, first, last+1) poly[m++] = pp[i];
    return m;
}


int main()
{
    while(scanf("%d", &n), n)
    {
        REP(i, n) scanf("%lf%lf", &p[i].x, &p[i].y);
        REP(i, n)
        {
            v[i] = p[(i+1)%n]-p[i];
            v2[i] = Normal(v[i]);
        }
        double l=0, r=20000, mid;
        while(r - l > eps)
        {
            mid = (l+r) / 2.0;
            REP(i, n) L[i] = Line(p[i]+v2[i]*mid, v[i]);
            int m = HalfplaneIntersection(L, n, poly);
            if(!m) r=mid; else l=mid;
        }
        printf("%.6f\n", l);
    }
    return 0;
}