POJ 3525 Most Distant Point from the Sea (半平面交求多边形中可以放入最大的圆的半径)

传送门：POJ 3525

题意：求多边形内到边上距离最远的距离。可以转化为求多边形内可以放入的最大的圆的半径。

前置技能： 半平面交求多边形中可以放入的最大的圆的半径。关于半平面交的知识可以参考我的博客：半平面交模板原理及应用。

思路：可以放入的最大圆的圆心（不是圆本身）一定在多边形的核内，别问我为什么，鬼才知道……我们可以假设，最大圆半径为 r 。如果我们可以将多边形的每条边向内推进 r 长度的距离，这时候可以想见多边形里面放不开这个圆了，并且圆心所在的多边形的核也一定是变成了一个点。否则就多边形的边就还可以向内推进，这个圆就不是最大的了。

当我们不知道最大圆半径的时候，我们可以用二分的方法求半径的值，如果对于当前假设的圆半径 r ，每条边向内推进 r 距离后，发现多边形的核不存在了，则说明当前圆半径太大了，反之，就是圆半径小了，用二分找到最大的圆半径即可。

代码：

给出两种代码，第二种是kuangbin模板上的，适用度更高，可以用作比赛时的模板。

//点是逆时针顺序给出的
//语言选择 c++  
/*
我们平移凸多边形的各个边，直到核不存在了，
平移的距离就是最大内切圆的半径。
*/
#include<stdio.h>
#include<string.h>
#include<iostream>
#include<algorithm>
#include<math.h>
#define eps 1e-10
#define MAXN 110
using namespace std;

int n,tol,bot,top; //bot和top指向双端队列的底部和顶部 
int dq[MAXN],order[MAXN]; //双端队列和存储排序后点的下标 

struct point
{ //点 
	double x,y;
} p[MAXN];

struct vct
{ //向量（直线）
	point start,end;
	double angle; //极角 
} v[MAXN],tmp[MAXN];

int dbcmp(double k)
{ //判断为0，还是为正负 
	if(fabs(k)<eps) return 0;
	return k>0? 1:-1;
}

double Multi(point p0,point p1,point p2)
{ //向量 p0p1 和 p0p2 相乘 
	return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);
}

int cmp(int x,int y)
{ //排序的比较函数 
	int d=dbcmp(v[x].angle-v[y].angle);
	if(!d) return dbcmp(Multi(v[x].start,v[y].start,v[y].end))>0;
	 //大于0取向量左半部分为半平面，小于0，取右半部分  
	return d<0;
}

void addline(double x1,double y1,double x2,double y2)
{ //增加一个向量（直线） 
	v[tol].start.x=x1;
	v[tol].start.y=y1;
	v[tol].end.x=x2;
	v[tol].end.y=y2;
	//atan2是求原点到(x,y)的直线与x轴的夹角，比atan稳定 
	v[tol].angle=atan2(y2-y1,x2-x1);
	//order[tol]=tol;
	tol++;
}

void GetIntersect(vct v1,vct v2,point &p)
{ //获取交点并返回给 p 
	double dot1,dot2;
	dot1=Multi(v1.start,v2.end,v1.end);
	dot2=Multi(v1.start,v2.start,v1.end);
	p.x=(dot1*v2.start.x-dot2*v2.end.x)/(dot1-dot2);		
	p.y=(dot1*v2.start.y-dot2*v2.end.y)/(dot1-dot2);
}

int Judge(vct v0,vct v1,vct v2)
{ //判断 v1 和 v2 的交点 pt 在半平面 v0内 
	point pt;
	GetIntersect(v1,v2,pt);
	return dbcmp(Multi(pt,v0.start,v0.end))<0;
	//大于0 pt在 v0 的左面，小于0 pt 在 v0 的右面
	//如果 pt 在半平面 v0 内则返回真 
}

int HalfPlaneIntersection()
{ //求解半平面交 
	int i,j;
	for(i=0;i<tol;i++) order[i]=i;
	sort(order,order+n,cmp); //极角排序 
	for(i=1,j=0;i<tol;i++) //排除极角相同的 
		if(dbcmp(tmp[order[i]].angle-tmp[order[j]].angle)>0)
			order[++j]=order[i]; 
	n=j+1; //个数 
	dq[0]=order[0]; //开始时入队两条直线 
	dq[1]=order[1];
	bot=0;
	top=1;
	for(i=2;i<tol;i++)
	{ //如果双端队列底端或顶端两个半平面的交点在当前半平面之外，则删除 
		while(bot<top&&Judge(tmp[order[i]],tmp[dq[top-1]],tmp[dq[top]]))
			top--;
		while(bot<top&&Judge(tmp[order[i]],tmp[dq[bot+1]],tmp[dq[bot]]))
			bot++;
		dq[++top]=order[i]; //将当前半平面加入队列的顶端 
	}
	while(bot<top&&Judge(tmp[dq[bot]],tmp[dq[top-1]],tmp[dq[top]])) top--;
	while(bot<top&&Judge(tmp[dq[top]],tmp[dq[bot+1]],tmp[dq[bot]])) bot++;	
	if(top-bot<=1) return 0;
	else return 1;
}

double GetDis(point a,point b)
{
	return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

void ChangePolygon(double h)
{ //将所有边往内缩进 h距离 
	int i;
	double dx,dy,len;
	for(i=0;i<tol;i++)
	{
		len=GetDis(v[i].start,v[i].end);
		dx=(v[i].start.y-v[i].end.y)/len*h;
		dy=(v[i].end.x-v[i].start.x)/len*h;
		tmp[i].start.x=v[i].start.x+dx;
		tmp[i].start.y=v[i].start.y+dy;
		tmp[i].end.x=v[i].end.x+dx;
		tmp[i].end.y=v[i].end.y+dy;
		tmp[i].angle=v[i].angle;
	}
}

double BSearch()
{ //二分搜索 
	double l=0,r=20000,mid;
	while(l+eps<r)
	{
		mid=(l+r)/2;
		ChangePolygon(mid); //将所有边往内缩进 mid距离 
		if(HalfPlaneIntersection()) l=mid; //如果存在半平面交
		else r=mid;
	}
	return l;
}

int main()
{
	int i;
	while(~scanf("%d",&n)&&n)
	{
		for(i=0;i<n;i++) scanf("%lf%lf",&p[i].x,&p[i].y);
		tol=0;
		for(i=0;i<n-1;i++) addline(p[i].x,p[i].y,p[i+1].x,p[i+1].y);		
		addline(p[i].x,p[i].y,p[0].x,p[0].y);
		printf("%.6lf\n",BSearch());
	}
	return 0;
}

#include<stdio.h>
#include<string.h>
#include<iostream>
#include<algorithm>
#include<math.h>
#define eps 1e-8
#define PI acos(-1.0)
#define MAXN 110
using namespace std;

int n;

int sgn(double x)
{ //符号函数
	if(fabs(x) < eps) return 0;
	if(x < 0) return -1;
	else return 1;
}

struct Point
{ //点 
	double x,y;
	Point(){}
	Point(double _x,double _y)
	{
		x = _x; y = _y;
	}
	Point operator -(const Point &b)const
	{
		return Point(x - b.x, y - b.y);
	}
	double operator ^(const Point &b)const
	{ //叉积
		return x*b.y - y*b.x;
	}
	double operator *(const Point &b)const
	{ //点积
		return x*b.x + y*b.y;
	}
};

struct Line
{ //向量 
	Point s,e; //两点 
	double k; //斜率 
	Line(){}
	Line(Point _s,Point _e)
	{ //构造
		s = _s; e = _e;
		k = atan2(e.y - s.y,e.x - s.x);
	}
	Point operator &(const Line &b)const
	{ //求两直线交点
		Point res = s;
		double t = ((s - b.s)^(b.s - b.e))/((s - e)^(b.s - b.e));
		res.x += (e.x - s.x)*t;
		res.y += (e.y - s.y)*t;
		return res;
	}
};

Line Q[MAXN];
Point p[MAXN]; //记录最初给的点集
Line line[MAXN]; //由最初的点集生成直线的集合
Point pp[MAXN]; //记录半平面交的结果的点集

//半平面交，直线的左边代表有效区域
bool HPIcmp(Line a,Line b)
{ //直线排序函数
	if(fabs(a.k - b.k) > eps)return a.k < b.k; //斜率排序
	//斜率相同我也不知道怎么办
	return ((a.s - b.s)^(b.e - b.s)) < 0;
}

void HPI(Line line[], int n, Point res[], int &resn)
{ //line是半平面交的直线的集合 n是直线的条数 res是结果
//的点集 resn是点集里面点的个数
	int tot = n;
	sort(line,line+n,HPIcmp);
	tot = 1;
	for(int i = 1;i < n;i++)
		if(fabs(line[i].k - line[i-1].k) > eps) //去掉斜率重复的
			line[tot++] = line[i];
	int head = 0, tail = 1;
	Q[0] = line[0];
	Q[1] = line[1];
	resn = 0;
	for(int i = 2; i < tot; i++)
	{
		if(fabs((Q[tail].e-Q[tail].s)^(Q[tail-1].e-Q[tail-1].s)) < eps ||
		fabs((Q[head].e-Q[head].s)^(Q[head+1].e-Q[head+1].s)) < eps)
			return;
		while(head < tail && (((Q[tail]&Q[tail-1]) -
		line[i].s)^(line[i].e-line[i].s)) > eps)
			tail--;
		while(head < tail && (((Q[head]&Q[head+1]) -
		line[i].s)^(line[i].e-line[i].s)) > eps)
			head++;
		Q[++tail] = line[i];
	}
	while(head < tail && (((Q[tail]&Q[tail-1]) -
	Q[head].s)^(Q[head].e-Q[head].s)) > eps)
		tail--;
	while(head < tail && (((Q[head]&Q[head-1]) -
	Q[tail].s)^(Q[tail].e-Q[tail].e)) > eps)
		head++;
	if(tail <= head + 1) return;
	for(int i = head; i < tail; i++)
		res[resn++] = Q[i]&Q[i+1];
	if(head < tail - 1)
		res[resn++] = Q[head]&Q[tail];
}

double dist(Point a,Point b)
{ //两点间距离
    return sqrt((a-b)*(a-b));
}

void change(Point a,Point b,Point &c,Point &d,double p)
{ //将线段ab往左移动距离p,修改得到线段cd
    double len=dist(a,b);
    /*三角形相似推出下面公式*/
    double dx=(a.y-b.y)*p/len;
    double dy=(b.x-a.x)*p/len;
    c.x=a.x+dx; c.y=a.y+dy;
    d.x=b.x+dx; d.y=b.y+dy;
}

double BSearch()
{ //二分搜索 
	double l=0,r=100000;
    double ans=0;
    while(r-l>=eps)
    {
        double mid=(l+r)/2;
        for(int i=0;i < n;i++)
        {
            Point t1,t2;
            change(p[i],p[(i+1)%n],t1,t2,mid);
            line[i]=Line(t1,t2);
        }
        int resn;
        HPI(line,n,pp,resn);
        //等于0说明移多了
        if(resn==0) r=mid-eps;
        else l=mid+eps;        
    }
    return l;
}

int main()
{
	while(~scanf("%d",&n)&&n)
	{		
		for(int i = 0;i < n;i++)
			scanf("%lf%lf",&p[i].x,&p[i].y);
		printf("%.6f\n",BSearch());
	}
	return 0;
}