AOJ 2121 : Castle Wall

問題リンク : Castle Wall | Aizu Online Judge

問題概要 :
単純多角形が与えられる
自分はこの多角形の2頂点を選んで新たに辺を付け加えることができる
これを0回以上行うことができるのだが、その際に追加した辺の長さの総和がr以下でなければならない
この処理が終わった後にできる多角形の面積を最大化し、その面積を出力せよ

解法 :

幾何 + 動的計画法
簡単な流れは、
1. 任意の２頂点を結んだ時にかかる辺の長さ、それによって得られる面積を求めグラフにしておく
2. それを使ってDP
となる

まず最初に与えられた多角形の頂点番号０番のものを凸包した際に残る頂点のいずれかにしておく
これは実際に凸包して探しても良いし、単純にループを回して多角形の点のうちもっとも左下にあるものを探しても良い
DPで頂点またぐ処理とか深く考えなくて良くなるのでやっておくと良い
凸包で残らない、つまり凹んだ部分にある点からDPを行うと０番を越えたりして面倒
凸包で残った頂点は別の点から辺を張る際にまたがれることがない

次に多角形に凸包をかけて凸多角形を得る
これらの凸多角形の頂点を順番に見ていき、もとの多角形の頂点番号で見て１つ以上飛んでいる(凹みがある）場所をみつける
見つけたら、もとの多角形の凹みとなっているその部分の中にある2点i,j( i != j && i < j, ただし j == 0 のとき i > j でも良い ) に辺を張った場合の情報をグラフに挿入しておく
この際に追加する線分が多角形に内包されていないか、別の辺と交差していないかを調べ、そのような場合は辺の追加をやめる

最後にDPで答えを求める
dp[頂点番号][面積] := 頂点番号までで面積をつくる際に必要な辺の長さの最小値
ただし、面積は本来の面積*2とする
今回全ての座標は整数なので、そのような多角形の面積を2倍すると整数になり扱いやすい
dp[i][j]がr以下のとき、jが答えの候補となるのでそのようなjのうち最大のものを出力すれば良い

コード :

#include<bits/stdc++.h>

#define REP(i,s,n) for(int i=s;i<n;i++)
#define rep(i,n) REP(i,0,n)
#define EPS (1e-8)
#define equals(a,b) (fabs((a)-(b)) < EPS)
#define COUNTER_CLOCKWISE 1
#define CLOCKWISE -1 
#define ONLINE_BACK 2
#define ONLINE_FRONT -2
#define ON_SEGMENT 0

using namespace std;

class Point{
public:
  double x,y;

  Point(double x = 0,double y = 0): x(x),y(y){}

  Point operator + (Point p){return Point(x+p.x,y+p.y);}
  Point operator - (Point p){return Point(x-p.x,y-p.y);}
  Point operator * (double a){return Point(a*x,a*y);}
  Point operator / (double a){return Point(x/a,y/a);}

  bool operator < (const Point& p) const{ return !equals(x,p.x)?x<p.x:(!equals(y,p.y)&&y<p.y); }
  bool operator == (const Point& p)const{ return fabs(x-p.x) < EPS && fabs(y-p.y) < EPS; }
};

struct Segment{
  Point p1,p2;
  Segment(Point p1 = Point(),Point p2 = Point()):p1(p1),p2(p2){}
  bool operator <  (const Segment& s) const { return ( p2 == s.p2 ) ? p1 < s.p1 : p2 < s.p2; }
  bool operator == (const Segment& s) const { return ( s.p1 == p1 && s.p2 == p2 ) || ( s.p1 == p2 && s.p2 == p1 ); }

};

typedef Point Vector;
typedef Segment Line;
typedef vector<Point> Polygon;

double dot(Point a,Point b){ return a.x*b.x + a.y*b.y; }
double cross(Point a,Point b){ return a.x*b.y - a.y*b.x; }
double norm(Point a){ return a.x*a.x+a.y*a.y; }
double abs(Point a){ return sqrt(norm(a)); }

int ccw(Point p0,Point p1,Point p2){
  Point a = p1-p0;
  Point b = p2-p0;
  if(cross(a,b) > EPS)return COUNTER_CLOCKWISE;
  if(cross(a,b) < -EPS)return CLOCKWISE;
  if(dot(a,b) < -EPS)return ONLINE_BACK;
  if(norm(a) < norm(b))return ONLINE_FRONT;
  return ON_SEGMENT;
}

bool intersectSS(Line s, Line t) {
  return ccw(s.p1,s.p2,t.p1)*ccw(s.p1,s.p2,t.p2) <= 0 &&
         ccw(t.p1,t.p2,s.p1)*ccw(t.p1,t.p2,s.p2) <= 0;
}

Point crosspoint(Line l, Line m) {
  double A = cross(l.p2 - l.p1, m.p2 - m.p1);
  double B = cross(l.p2 - l.p1, l.p2 - m.p1);
  if (abs(A) < EPS && abs(B) < EPS) return m.p1; // same line
  if (abs(A) < EPS) assert(false); // !!!PRECONDITION NOT SATISFIED!!!
  return m.p1 +   (m.p2 - m.p1) * (B / A);
}

double cross3p(Point p,Point q,Point r) { return (r.x-q.x) * (p.y -q.y) - (r.y - q.y) * (p.x - q.x); }
bool collinear(Point p,Point q,Point r) { return fabs(cross3p(p,q,r)) < EPS; }
bool onSegment(Point p,Point q,Point r){
  return collinear(p,q,r) && equals(sqrt(pow(p.x-r.x,2)+pow(p.y-r.y,2)) + sqrt(pow(r.x-q.x,2) + pow(r.y-q.y,2) ),sqrt(pow(p.x-q.x,2)+pow(p.y-q.y,2)) ) ;
}

double angle(Point a,Point b,Point c) {
  double ux = b.x - a.x, uy = b.y - a.y;
  double vx = c.x - a.x, vy = c.y - a.y;
  return acos((ux*vx + uy*vy)/sqrt((ux*ux + uy*uy) * (vx*vx + vy*vy)));
}  
  
bool inPolygon(Polygon poly,Point p){
  if((int)poly.size() == 0)return false;
  rep(i,poly.size()) if(onSegment(poly[i],poly[(i+1)%poly.size()],p))return true;
  double sum = 0;
  for(int i=0; i < (int)poly.size() ;i++) {
    if( equals(cross(poly[i]-p,poly[(i+1)%poly.size()]-p),0.0) ) continue;
    if( cross3p(p,poly[i],poly[(i+1)%poly.size()]) < 0 ) sum -= angle(p,poly[i],poly[(i+1)%poly.size()]);
    else                                                 sum += angle(p,poly[i],poly[(i+1)%poly.size()]);
  }
  const double eps = 1e-6;
  return (fabs(sum - 2*M_PI) < eps || fabs(sum + 2*M_PI) < eps);
}  

Polygon andrewScan(Polygon s){
  Polygon u,l;
  if(s.size() < 3)return s;
  sort(s.begin(),s.end());
  u.push_back(s[0]); u.push_back(s[1]);
  l.push_back(s[s.size()-1]); l.push_back(s[s.size()-2]);

  for(int i=2;i<s.size();i++) {
    for(int n=u.size();n >= 2 && ccw(u[n-2],u[n-1],s[i]) != CLOCKWISE; n--) u.pop_back();
    u.push_back(s[i]);
  }

  for(int i=s.size()-3; i>=0 ; i--){
      for(int n=l.size(); n >= 2 && ccw(l[n-2],l[n-1],s[i]) != CLOCKWISE; n--) l.pop_back();
      l.push_back(s[i]);
  }
  reverse(l.begin(),l.end());
  for(int i = u.size()-2; i >= 1; i--) l.push_back(u[i]);
  return l;
}

bool inPolygon(const Polygon &poly,Segment seg){
  vector<Point> vp;
  vp.push_back(seg.p1);
  vp.push_back(seg.p2);
  rep(i,(int)poly.size()) {
    Segment edge = Segment(poly[i],poly[(i+1)%(int)poly.size()]);
    if( equals(cross(seg.p1-seg.p2,edge.p1-edge.p2),0.0) ) {
      if( onSegment(seg.p1,seg.p2,edge.p1) ) vp.push_back(edge.p1);
      if( onSegment(seg.p1,seg.p2,edge.p2) ) vp.push_back(edge.p2);
    } else {
      if( intersectSS(seg,edge) ) vp.push_back(crosspoint(seg,edge));
    }
  }
  sort(vp.begin(),vp.end());
  vp.erase(unique(vp.begin(),vp.end()),vp.end());
  rep(i,(int)vp.size()-1) {
    Point middle_point = ( vp[i] + vp[i+1] ) / 2.0;
    if( !inPolygon(poly,middle_point) ) return false;
  }
  return true;  
}

double getArea(const Polygon &poly){
  double sum = 0;
  rep(i,(int)poly.size()) sum += cross(poly[i],poly[(i+1)%(int)poly.size()]);
  return fabs(sum) * 0.5;
}

struct Edge {
  int dst,area;
  double length;
};

const double DINF = 1e30;

bool LTE(double a,double b) { return equals(a,b) || a < b; } 
bool LT(double a,double b) { return !equals(a,b) && a < b; }

bool veri(const Polygon& poly,Segment seg) {
  rep(i,(int)poly.size()) {
    Segment edge = Segment(poly[i],poly[(i+1)%poly.size()]);
    if( edge == seg ) return false;
    vector<Point> vec;
    vec.push_back(seg.p1), vec.push_back(seg.p2), vec.push_back(edge.p1), vec.push_back(edge.p2);
    sort(vec.begin(),vec.end());
    vec.erase(unique(vec.begin(),vec.end()),vec.end());
    if( vec.size() == 3 ) continue;
    assert( vec.size() == 4);
    if( equals(cross(edge.p1-edge.p2,seg.p1-seg.p2),0.0) ) {
      if( onSegment(edge.p1,edge.p2,seg.p1) || onSegment(edge.p1,edge.p2,seg.p2) ||
          onSegment(seg.p1,seg.p2,edge.p1)  || onSegment(seg.p1,seg.p2,edge.p2) ) return false;
    } else {
      if( intersectSS(edge,seg) ) return false;
    }
  }
  return true;
}

int V;
Polygon poly;
double r;
const int LIMIT = 21000;
double dp[70][LIMIT+1];

vector<Edge> G[70];

void make_graph(){
  rep(i,poly.size()+1) G[i].clear();
  Polygon convex = andrewScan(poly);
  map<Point,int> mp;
  rep(i,(int)poly.size()) mp[poly[i]] = i;
  rep(i,(int)convex.size()) {
    int s = mp[convex[i]], t = mp[convex[(i+1)%(int)convex.size()]];
    if( (s+1)%V == t ) continue;
    for(int j=s;;(j+=1)%=V){
      for(int l=j;;(l+=1)%=V){
        if( j != l ){
          Segment seg = Segment(poly[j],poly[l]);
          if( !inPolygon(poly,seg) && veri(poly,seg) ) {
            vector<Point> vec;
            for(int k=j;;(k+=1)%=V) {
              vec.push_back(poly[k]);
              if( k == l ) break;
            }
            if( l ) G[min(j,l)].push_back((Edge){max(l,j),getArea(vec)*2,abs(poly[j]-poly[l])});
            else  G[j].push_back((Edge){V,getArea(vec)*2,abs(poly[j]-poly[l])});
          }
        }
        if( l == t ) break;
      }
      if( j == t ) break;
    }
  }
}

int CNT;

void compute(){
  make_graph();
  rep(i,V+1) rep(j,LIMIT+1) dp[i][j] = DINF;
  double maxi = getArea(poly);
  dp[0][(int)(maxi*2)] = 0;
  rep(i,V){
    rep(area,LIMIT)if( !equals(dp[i][area],DINF) ){
      rep(j,(int)G[i].size()){
        int next = G[i][j].dst;
        int cost = G[i][j].area;
        double length = G[i][j].length;
        if( area + cost < LIMIT && LT(dp[i][area]+length,dp[next][area+cost]) ) {
          dp[next][area+cost] = dp[i][area]+length;
        }
      }
    }
    rep(j,LIMIT) dp[i+1][j] = min(dp[i+1][j],dp[i][j]);
  }
  for(int i=LIMIT;i>=0;i--) if( !equals(dp[V][i],DINF) && LTE(dp[V][i],r) ) { 
      maxi = i * 0.5; break; 
    }
  printf("case %d: %.1f\n",CNT++,maxi);
}

int main(){
  int temp;
  CNT = 1;
  while( cin >> V >> temp, V|temp){
    r = temp;
    poly.clear();
    poly.resize(V);
    rep(i,V) cin >> poly[i].x >> poly[i].y;
    Point mini = poly[0];
    rep(i,V) mini = min(mini,poly[i]);
    int sp = -1;
    rep(i,V) if( mini == poly[i] ) { sp = i; break; }
    Polygon ps;
    for(int i=sp;i<V;i++) ps.push_back(poly[i]);
    for(int i=0;i<sp;i++) ps.push_back(poly[i]);
    poly = ps; 
    compute();
  }
  return 0;
}