poj 3525 Most Distant Point from the Sea - 求到海岸最远的点 - 半平面交

/*
poj 3525 Most Distant Point from the Sea - 求到海岸最远的点 - 半平面交

就是求多边形最大的内接圆的半径

枚举半径
将所有海岸沿法向量向里推进半径的距离，求按平面的交，若多边形还有核，则这个半径存在
*/
#include<stdio.h>
#include<math.h>
#include <algorithm>
using namespace std;

const double eps=1e-9;
struct point
{
double x,y;
point(){}
point(double a,double b):x(a),y(b){}
}dian[803];
point jiao[803];
struct line
{
point s,e;
double angle;
}xian[803];
int n,yong;
bool mo_ee(double x,double y)
{
double ret=x-y;
if(ret<0) ret=-ret;
if(ret<eps) return 1;
return 0;
}
bool mo_gg(double x,double y)  {   return x > y + eps;} // x > y
bool mo_ll(double x,double y)  {   return x < y - eps;} // x < y
bool mo_ge(double x,double y) {   return x > y - eps;} // x >= y
bool mo_le(double x,double y) {   return x < y + eps;}     // x <= y
point mo_intersection(point u1,point u2,point v1,point v2)
{
point ret=u1;
double t=((u1.x-v1.x)*(v1.y-v2.y)-(u1.y-v1.y)*(v1.x-v2.x))
/((u1.x-u2.x)*(v1.y-v2.y)-(u1.y-u2.y)*(v1.x-v2.x));
ret.x+=(u2.x-u1.x)*t;
ret.y+=(u2.y-u1.y)*t;
return ret;
}
double mo_xmult(point p2,point p0,point p1)//p1在p2左返回负，在右边返回正
{
return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);
}

void mo_HPI_addl(point a,point b)
{
xian[yong].s=a;
xian[yong].e=b;
xian[yong].angle=atan2(b.y-a.y,b.x-a.x);
yong++;
}
//半平面交
bool mo_HPI_cmp(const line& a,const line& b)
{
if(mo_ee(a.angle,b.angle))
{
return mo_gg( mo_xmult(b.e,a.s,b.s),0);
}else
{
return mo_ll(a.angle,b.angle);
}
}
int mo_HPI_dq[803];
bool mo_HPI_isout(line cur,line top,line top_1)
{
point jiao=mo_intersection(top.s,top.e,top_1.s,top_1.e);
return mo_ll( mo_xmult(cur.e,jiao,cur.s),0);//若顺时针时应为mo_gg
}
int mo_HalfPlaneIntersect(line *xian,int n,point *jiao)
{
int i,j,ret=0;
sort(xian,xian+n,mo_HPI_cmp);
for (i = 0, j = 0; i < n; i++)
{
if (mo_gg(xian[i].angle,xian[j].angle))
{
xian[++j] = xian[i];
}
}
n=j+1;
mo_HPI_dq[0]=0;
mo_HPI_dq[1]=1;
int top=1,bot=0;
for (i = 2; i < n; i++)
{
while (top > bot && mo_HPI_isout(xian[i], xian[mo_HPI_dq[top]], xian[mo_HPI_dq[top-1]])) top--;
while (top > bot && mo_HPI_isout(xian[i], xian[mo_HPI_dq[bot]], xian[mo_HPI_dq[bot+1]])) bot++;
mo_HPI_dq[++top] = i; //当前半平面入栈
}
while (top > bot && mo_HPI_isout(xian[mo_HPI_dq[bot]], xian[mo_HPI_dq[top]], xian[mo_HPI_dq[top-1]])) top--;
while (top > bot && mo_HPI_isout(xian[mo_HPI_dq[top]], xian[mo_HPI_dq[bot]], xian[mo_HPI_dq[bot+1]])) bot++;
mo_HPI_dq[++top] = mo_HPI_dq[bot];
for (ret = 0, i = bot; i < top; i++, ret++)
{
jiao[ret]=mo_intersection(xian[mo_HPI_dq[i+1]].s,xian[mo_HPI_dq[i+1]].e,xian[mo_HPI_dq[i]].s,xian[mo_HPI_dq[i]].e);
}
return ret;
}
//求法向量
point mo_getfaxian(point xiang)
{
point a;
if(mo_ee(xiang.x,0))
{
a.x=1;
a.y=0;
return a;
}else if(mo_ee(xiang.y,0))
{
a.x=0;
a.y=1;
return a;
}else
{
a.x=1;
a.y=-1.0*xiang.x/xiang.y;
return a;
}

}
//线se推进一段距离
line mo_line_tui(point s,point e,double tui)
{
line ret;
point fa=mo_getfaxian(point(e.x-s.x,e.y-s.y));
if(mo_ll(mo_xmult(e,point(s.x+fa.x,s.y+fa.y),s),0))
{
fa.x=-fa.x;
fa.y=-fa.y;
}
double len=sqrt(fa.x*fa.x+fa.y*fa.y);
len=len/tui;
fa.x=fa.x/len;
fa.y=fa.y/len;

s.x=s.x+fa.x,s.y=s.y+fa.y;
e.x=e.x+fa.x,e.y=e.y+fa.y;

ret.s=s;
ret.e=e;
return ret;
}
int kexing(double jia)
{
int i;
yong=0;
line temp;
for(i=0;i<n;++i)
{
temp=mo_line_tui(dian[i],dian[(i+1)%n],jia);
mo_HPI_addl(temp.s,temp.e);
}
int ret=mo_HalfPlaneIntersect(xian,n,jiao);
if(ret<3)
{
return 0;
}
return 1;
}
int main()
{
int i;
double banjing;
while(scanf("%d",&n),n)
{
yong=0;
for(i=0;i<n;++i)
{
scanf("%lf%lf",&dian[i].x,&dian[i].y);
}
double shang=5001,xia=0,mid;
while(xia+1e-6<shang)
{
mid=(shang+xia)/2;
int ret=kexing(mid);
if(ret)
{
banjing=mid;
xia=mid+1e-6;
}else
{
shang=mid-1e-6;
}
}
printf("%lf\n",banjing);
}
return 0;
}