Most Distant Point from the Sea UVA - 1396

半平面交模板题

#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int N = 1e6+10;
const double eps = 1e-10;//精度
//atan2(y,x)求极角的弧度
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)//判断位置，先按照x坐标排序然后按照y坐标排序
{
    return a.x < b.x || (a.x == b.x && a.y < b.y);
}

int dcmp(double x)//判断x的符号
{
    if(fabs(x) < eps)
        return 0;
    else
        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));
}
//A×B,B在A的左边是正，在右边是负
double Cross(Vector A,Vector B)//两个向量的叉积
{
    return A.x*B.y - A.y*B.x;
}
//把向量向左旋转90°,得到法向量,得保证A向量不是零向量
Vector Normal(Vector A)
{
    double L = Length(A);
    return Vector(-A.y/L,A.x/L);
}
//两条直线的交点,向量式--参数方程.前提是v,w向量不平行,Cross(v,w)!=0
struct Line{
    Point P;//直线上任意一点
    Vector v;//方向向量，他的左边对应半平面
    double ang;//极角，即从x正半轴旋转到向量v所需要的角
    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在有向直线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;

}
//求半平面交，答案放入poly,return 交成多边形的交点数目
int HalfplaneIntersection(Line *L,int n,Point *poly)
{
    sort(L,L+n);//极角排序
    int first,last;//双端队列的第一个元素和最后一个元素的下标
    Point *p = new Point[n];//p[i]为q[i]和q[i+1]的下标
    Line *q = new Line[n];//双端队列
    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;
}
void solve()
{
    Point p[200],poly[200];
    Line L[200];
    Vector v[200],v2[200];//v2的是单位法向量
    int n;
    while(cin >> n && n)
    {
        int m,x,y;
        for(int i=0;i<n;i++)
        {
            cin >> x >> y;
            p[i]= Point(x,y);
        }
        for(int i=0;i<n;i++)
        {
            v[i] = p[(i+1)%n] - p[i];
            v2[i] = Normal(v[i]);//求出单位法向量
        }
        double left = 0,right = 20000;
        while(right - left > 1e-7)
        {
            double mid = left + (right -left) / 2;
            for(int i=0;i<n;i++)
            {
                L[i] = Line(p[i]+v2[i]*mid,v[i]);//往法向量方向移动mid个单位,收缩多边形
            }
            m = HalfplaneIntersection(L,n,poly);
            if(!m)
                right = mid;
            else
                left = mid;
        }
        printf("%0.6lf\n",left);
    }
}
int main()
{
    int T = 1;
//    cin >> T;
    while(T--)
    {
        solve();
    }
    return 0;
}