  code
#include<iostream>
#include<cmath>
using namespace std;
const double eps=1e-10;
const double limit=1e-10;
struct PT
{
double x,y;
};
struct LN
{
double a,b,c;
};
struct PL
{
int n;
PT p[1501];
};
PL poly,newpoly,tmpoly;
int cas;
int n;
LN getline(PT a,PT b)
{
LN l;
l.a=b.y-a.y;
l.b=a.x-b.x;
l.c=b.x*a.y-a.x*b.y;
return l;
}
PT getinterp(LN l1,LN l2)
{
PT ans;
if(fabs(l1.b)<eps)
{
ans.x=-l1.c/l1.a;
ans.y=(-l2.c-l2.a*ans.x)/l2.b;
}
else
{
ans.x=(l2.c*l1.b-l1.c*l2.b)/(l1.a*l2.b-l2.a*l1.b);
ans.y=(-l1.c-l1.a*ans.x)/l1.b;
}
return ans;
}
double Cha(PT a,PT b,PT c)
{
return (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x);
}
PT change(PT s,PT e,PT st,double L)
{
PT tmp;
tmp.x=-(e.y-s.y);
tmp.y=e.x-s.x;
double len=sqrt((e.x-s.x)*(e.x-s.x)+(e.y-s.y)*(e.y-s.y));
tmp.x/=len;tmp.x*=L;
tmp.y/=len;tmp.y*=L;
tmp.x+=st.x;
tmp.y+=st.y;
return tmp;
}
bool jud(double h)
{
int i,j;
newpoly=poly;
for(i=0;i<poly.n;i++)
{
tmpoly.n=0;
PT s,e;
s=change(poly.p[i],poly.p[(i+1)%poly.n],poly.p[i],h);
e=change(poly.p[i],poly.p[(i+1)%poly.n],poly.p[(i+1)%poly.n],h);
for(j=0;j<newpoly.n;j++)
{
double r1,r2;
r1=Cha(s,e,newpoly.p[j]);
r2=Cha(s,e,newpoly.p[(j+1)%newpoly.n]);
if(r1>=0)
{
tmpoly.p[tmpoly.n++]=newpoly.p[j];
}
if(r1*r2<0)
{
LN l1=getline(s,e);
LN l2=getline(newpoly.p[j],newpoly.p[(j+1)%newpoly.n]);
PT interp=getinterp(l1,l2);
tmpoly.p[tmpoly.n++]=interp;
}
}
newpoly=tmpoly;
}
if(newpoly.n)
return 1;
return 0;
}
int main()
{
int i,j;
while(1)
{
scanf("%d",&poly.n);
if(!poly.n)
break;
for(i=0;i<poly.n;i++)
scanf("%lf%lf",&poly.p[i].x,&poly.p[i].y);
double left=0,right=100000000;
while(left+limit<right)
{
double mid=(left+right)/2;
if(jud(mid))
left=mid;
else right=mid;
}
printf("%.6lf\n",right);
}
}
求凸多边形最大内切圆半径。 方法就是把每条边向内推进R，对得到的新的边集进行半平面交，看是否得到空集。 R用二分枚举得到。 半平面交用的O(n^2)的。
|