AOJ 1283 Most Distant Point from the Sea

問題

http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1283
多角形の内部の点で、辺との距離の最小値が最大になる点を求めよ。

解法

初めは他の人と同じように想定解の二分探索でやっていたが、謎のWAが解決できなかったから次のような山登り法?のような探索に変えた。

適当な内点を初期値として、最も距離が小さい辺と逆側に少しずつ動かしていった。
このアイデアは http://d.hatena.ne.jp/y_mazun/20121218/1355851934 を参考にした。
適当な3点を選び、その内心を最初の点としている。そこから最も距離が近い辺を探し、その方向とは逆向きのベクトルを足す。足すことによって点が円の外に出てしまうならその長さを足すことを諦めて、足す長さを短くする。これを答えに求められる精度が出る回数だけ繰り返す。 SpaghettiSourceのライブラリを使えば多角形の内外判定がlogの時間でできるのでだいたい となるが、これにかかる定数で精度が決まるのでそこを調節する必要がある。

最初の点の選び方によって失敗することがあるようで、少しランダム要素を加えて3回に2回以上はACされるように調整した。

ソースコード

double const eps=1e-8;
double const inf = 1e10;
double const pi = acos(-1);

typedef complex<double> P;
inline double cross(P const&a, P const&b){return imag(conj(a)*b);}
inline double dot(P const&a, P const&b){return real(conj(a)*b);}
struct L : public vector<P> {
    L(const P &a, const P &b){push_back(a); push_back(b);}
    L(){}
};
P projection(L const&l, P const&p){double t = dot(p-l[0],l[0]-l[1])/norm(l[0]-l[1]);return l[0]+t*(l[0]-l[1]);}
inline double distanceLP(L const&l, P const&p){return abs(p-projection(l, p));}

P innerCenter(P const&a, P const&b, P const&c){
    double la=abs(b-c), lb=abs(c-a), lc=abs(a-b);
    return (la*a+lb*b+lc*c)/(la+lb+lc);
}

int n;
P ps[128];
L es[128];

enum { OUT, ON, IN };
int convex_contains(P const&p) {
    P g=(ps[0]+ps[n/3]+ps[2*n/3])/3.0;
    int a=0,b=n;
    while (a+1 < b) {
        int c=(a+b)/2;
        if (cross(ps[a]-g, ps[c]-g) > 0) {
            if (cross(ps[a]-g, p-g) > 0 && cross(ps[c]-g, p-g) < 0) b = c;
            else                                                    a = c;
        } else {
            if (cross(ps[a]-g, p-g) < 0 && cross(ps[c]-g, p-g) > 0) a = c;
            else                                                    b = c;
        }
    }
    b %= n;
    if (cross(ps[a] - p, ps[b] - p) < 0) return 0;
    if (cross(ps[a] - p, ps[b] - p) > 0) return 2;
    return 1;
}

double solve(){
    rotate(ps,ps+rand()%n,ps+n);
    ps[n]=ps[0];
    rep(i,n){
        es[i]={ps[i],ps[i+1]};
    }
    P p=innerCenter(ps[0],ps[n/3],ps[n/3*2]);

    if(n==3)return distanceLP(es[0],p);

    double const r[]={100,10,1,1e-1,1e-2,1e-3,1e-4,1e-5,1e-6,1e-7,1e-8};
    double static d[128];
    int near=0;
    rep(t,11){
        rep(u,3000){
            rep(i,n)d[i]=distanceLP(es[i],p);
            near=0;
            rep(i,n)
                if(d[i]<d[near]){near=i;
            }
            P h=projection(es[near],p);
            P np=p-(h-p)/d[near]*r[t];
            if(convex_contains(np)==IN){
                swap(p,np);
            } else {
                break;
            }
        }
    }
    return d[near];
}

int main(){
    srand(time(0));
    while(cin>>n && n){
        rep(i,n){
            double x,y;
            scanf("%lf%lf",&x,&y);
            ps[i]=P(x,y);
        }
        printf("%lf\n",solve());
    }
}

幾何の問題は完全にSpagettiSourceに頼り切っているので自力でライブラリ書かないといけないと思った。