392      阿裏雲  技術社區[雲棲]

POJ 3525 二分+半麵相交

題意：給定一個凸多邊形問能包裹圓的最大半徑。

將每條邊向內移動距離d，知道半麵相交的內核不存在，確定最大的d並且達到精度就可以了。

#include<cmath>
#include<cstdio>
#include<iostream>
#include<algorithm>
using namespace std;
const int mm=111;
typedef double DIY;
struct point
{
    DIY x,y;
    point() {}
    point(DIY _x,DIY _y):x(_x),y(_y) {}
} g[mm];
point MakeVector(point &P,point &Q)
{
    return point(Q.x-P.x,Q.y-P.y);
}
DIY CrossProduct(point P,point Q)
{
    return P.x*Q.y-P.y*Q.x;
}
DIY MultiCross(point P,point Q,point R)
{
    return CrossProduct(MakeVector(Q,P),MakeVector(Q,R));
}
struct halfPlane
{
    point s,t;
    double angle;
    halfPlane() {}
    halfPlane(point _s,point _t):s(_s),t(_t) {}
    halfPlane(DIY sx,DIY sy,DIY tx,DIY ty):s(sx,sy),t(tx,ty) {}
    void GetAngle()
    {
        angle=atan2(t.y-s.y,t.x-s.x);
    }
} hp[mm],q[mm];
point IntersectPoint(halfPlane P,halfPlane Q)
{
    DIY a1=CrossProduct(MakeVector(P.s,Q.t),MakeVector(P.s,Q.s));
    DIY a2=CrossProduct(MakeVector(P.t,Q.s),MakeVector(P.t,Q.t));
    return point((P.s.x*a2+P.t.x*a1)/(a2+a1),(P.s.y*a2+P.t.y*a1)/(a2+a1));
}
bool cmp(halfPlane P,halfPlane Q)
{
    if(fabs(P.angle-Q.angle)<1e-8)
        return MultiCross(P.s,P.t,Q.s)>0;
    return P.angle<Q.angle;
}
bool IsParallel(halfPlane P,halfPlane Q)
{
    return fabs(CrossProduct(MakeVector(P.s,P.t),MakeVector(Q.s,Q.t)))<1e-8;
}
void HalfPlaneIntersect(int n,int &m)
{
    sort(hp,hp+n,cmp);
    int i,l=0,r=1;
    for(m=i=1; i<n; ++i)
        if(hp[i].angle-hp[i-1].angle>1e-8)hp[m++]=hp[i];
    n=m;
    m=0;
    q[0]=hp[0],q[1]=hp[1];
    for(i=2; i<n; ++i)
    {
        if(IsParallel(q[r],q[r-1])||IsParallel(q[l],q[l+1]))return;
        while(l<r&&MultiCross(hp[i].s,hp[i].t,IntersectPoint(q[r],q[r-1]))>0)--r;
        while(l<r&&MultiCross(hp[i].s,hp[i].t,IntersectPoint(q[l],q[l+1]))>0)++l;
        q[++r]=hp[i];
    }
    while(l<r&&MultiCross(q[l].s,q[l].t,IntersectPoint(q[r],q[r-1]))>0)--r;
    while(l<r&&MultiCross(q[r].s,q[r].t,IntersectPoint(q[l],q[l+1]))>0)++l;
    q[++r]=q[l];
    for(i=l; i<r; ++i)
        g[m++]=IntersectPoint(q[i],q[i+1]);
}
point data[105];
int n,m;
void getFXXL(point a,point b,double &ax,double &ay)  //求垂直該向量的向右方向的方向向量
{
    double m=b.x-a.x,n=b.y-a.y;
    ax=n/sqrt(m*m+n*n),ay=(-m)/sqrt(m*m+n*n);
}
bool slove(double d)    //判斷是否存在內核
{
    double ax,ay;
    for(int i=0; i<n; i++)
    {
        getFXXL(data[(i+1)%n],data[i],ax,ay);
        g[i].x=data[i].x+ax*d,g[i].y=data[i].y+ay*d;
        g[(i+1)%n].x=data[(i+1)%n].x+ax*d,g[(i+1)%n].y=data[(i+1)%n].y+ay*d;
        hp[i]=halfPlane(g[i],g[(i+1)%n]);
        hp[i].GetAngle();
    }
    HalfPlaneIntersect(n,m);
    if(m>2)
        return 1;
    return 0;
}
int main()
{
    double d;
    while(~scanf("%d",&n),n)
    {
        for(int i=0; i<n; i++)
            scanf("%lf%lf",&data[i].x,&data[i].y);
        double l=0,r=1000000;
        while(r>l)         //二分求d
        {
            d=(r+l)/2;
            bool a1=slove(d),a2=slove(d+1e-8);
            if(a1&&!a2)
                break;
            else if(!a1)
                r=d;
            else
                l=d;
        }
        printf("%.6f\n",d);
    }
    return 0;
}

最後更新：2017-04-03 21:30:13

  上一篇： 機架服務器 遠程電源管理 備忘

  下一篇： 推薦一個IE6下js調試工具（Companion.JS）