AOJ 2154: Shore Erosion - SSRS の競プロ日記

この記事は帰ってきた AOJ-ICPC Advent Calendar 2022 22 日目の記事です。

https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2154

解法

与えられた多角形の内部を $D$ 、頂点を反時計回りに $P _ 0, P _ 1, \dots, P _ {N-1}$ とします。また、 $P _ N = P _ 0$ であるとします。

辺 $P _ i P _ {i+1}$ からのマンハッタン距離が $R$ 以下である点の集合を $D _ i$ とすると、島のある点が浸食される条件は、その点がある $i$ に対し $D _ i$ に含まれることです。よって、浸食された後の島は $D \setminus (D _ 0 \cup D _ 1 \cup D _ {n-1})$ となります。

領域 $D _ i$ は四角形もしくは六角形になります。具体的には、点 $A ( A _ x , A _ y ) , B ( B _ x , B _ y )$ があるとき、線分 $AB$ からのマンハッタン距離が $R$ 以下である領域は以下のように得られます。

$A _ x + A _ y \leq B _ x + B _ y$ と仮定しても一般性を失わないので、そのように仮定する。
$A _ x - A _ y > B _ x - B _ y$ の場合、 $6$ 点 $(A _ x - R, A _ y), (B _ x - R, B _ y), (B _ x, B _ y + R), (B _ x + R, B _ y), (A _ x + R, A _ y), (A _ x, A _ y - R)$ を頂点とする六角形である。
そうでない場合、 $6$ 点 $(A _ x - R, A _ y), (A _ x, A _ y + R), (B _ x, B _ y + R), (B _ x + R, B _ y), (B _ x, B _ y - R), (A _ x, A _ y - R)$ を頂点とする六角形である。

Slab Decomposition で平面を $D, D _ 0, D _ 1, \dots, D _ {N-1}$ の辺で分割し、それぞれの領域が浸食された後の島であるかどうかを調べます。ビット演算を用いて、それぞれの領域がどれに含まれるかを計算すると得られます。その後、島である領域と島でない領域が隣接する部分それぞれについて、長さを足し合わせればよいです。

実装

judge.u-aizu.ac.jp

#include <bits/stdc++.h>
using namespace std;
const double eps = 0.000001;
int sign(double x){
  if (x > eps){
    return 1;
  } else if (x < -eps){
    return -1;
  } else {
    return 0;
  }
}
struct point{
  double x, y;
  point(){
  }
  point(double x, double y): 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 k){
    return point(x * k, y * k);
  }
  point operator /(double k){
    return point(x / k, y / k);
  }
};
double abs(point P){
  return sqrt(P.x * P.x + P.y * P.y);
}
double dist(point P, point Q){
  return abs(Q - P);
}
point rotate(point P, double t){
  return point(P.x * cos(t) - P.y * sin(t), P.x * sin(t) + P.y * cos(t));
}
double cross(point P, point Q){
  return P.x * Q.y - P.y * Q.x;
}
struct line{
  point A, B;
  line(){
  }
  line(point A, point B): A(A), B(B){
  }
};
point vec(line L){
  return L.B - L.A;
}
bool is_parallel(line L1, line L2){
  return abs(cross(vec(L1), vec(L2))) < eps;
}
bool segment_line_intersect(line L1, line L2){
  return sign(cross(L1.A - L2.A, vec(L2))) * sign(cross(L1.B - L2.A, vec(L2))) <= 0;
}
bool segment_intersect(line L1, line L2){
  return segment_line_intersect(L1, L2) && segment_line_intersect(L2, L1);
}
point line_intersection(line L1, line L2){
  return L1.A + vec(L1) * cross(L2.A - L1.A, vec(L2)) / cross(vec(L1), vec(L2));
}
vector<double> unique_eps(vector<double> A){
  int N = A.size();
  sort(A.begin(), A.end());
  vector<double> A2;
  A2.push_back(A[0]);
  for (int i = 1; i < N; i++){
    if (A[i] - A[i - 1] > eps){
      A2.push_back(A[i]);
    }
  }
  return A2;
}
template <typename T>
struct slab_decomposition{
  function<T(T, T)> f;
  T e;
  vector<pair<line, T>> S;
  vector<double> X;
  int N;
  vector<vector<double>> Y;
  vector<int> P;
  vector<vector<tuple<int, int, T>>> E;
  slab_decomposition(function<T(T, T)> f, T e): f(f), e(e){
  }
  void add_segment(line s, T x){
    S.push_back(make_pair(s, x));
  }
  vector<vector<pair<T, int>>> build(){
    int M = S.size();
    for (int i = 0; i < M; i++){
      if (S[i].first.A.x > S[i].first.B.x){
        swap(S[i].first.A, S[i].first.B);
      }
    }
    for (int i = 0; i < M; i++){
      X.push_back(S[i].first.A.x);
      X.push_back(S[i].first.B.x);
    }
    for (int i = 0; i < M; i++){
      for (int j = i + 1; j < M; j++){
        if (!is_parallel(S[i].first, S[j].first) && segment_intersect(S[i].first, S[j].first)){
          X.push_back(line_intersection(S[i].first, S[j].first).x);
        }
      }
    }
    X = unique_eps(X);
    N = X.size();
    vector<int> L(M), R(M);
    for (int i = 0; i < M; i++){
      L[i] = lower_bound(X.begin(), X.end(), S[i].first.A.x - eps) - X.begin();
      R[i] = lower_bound(X.begin(), X.end(), S[i].first.B.x - eps) - X.begin();
    }
    vector<vector<double>> A(M, vector<double>(N));
    Y.resize(N);
    for (int i = 0; i < M; i++){
      for (int j = L[i]; j <= R[i]; j++){
        A[i][j] = line_intersection(S[i].first, line(point(X[j], 0), point(X[j], 1))).y;
        Y[j].push_back(A[i][j]);
      }
    }
    for (int i = 0; i < N; i++){
      Y[i] = unique_eps(Y[i]);
    }
    vector<vector<int>> B(M, vector<int>(N));
    for (int i = 0; i < M; i++){
      for (int j = L[i]; j <= R[i]; j++){
        B[i][j] = lower_bound(Y[j].begin(), Y[j].end(), A[i][j] - eps) - Y[j].begin();
      }
    }
    E.resize(N - 1);
    for (int i = 0; i < N - 1; i++){
      E[i].push_back(make_tuple(-1, -1, e));
    }
    for (int i = 0; i < M; i++){
      for (int j = L[i]; j < R[i]; j++){
        E[j].push_back(make_tuple(B[i][j], B[i][j + 1], S[i].second));
      }
    }
    for (int i = 0; i < N - 1; i++){
      E[i].push_back(make_tuple(Y[i].size(), Y[i + 1].size(), e));
    }
    for (int i = 0; i < N - 1; i++){
      sort(E[i].begin(), E[i].end());
      if (!E[i].empty()){
        vector<tuple<int, int, T>> E2;
        int cnt = E[i].size();
        E2.push_back(E[i][0]);
        for (int j = 1; j < cnt; j++){
          if (get<0>(E[i][j - 1]) == get<0>(E[i][j]) && get<1>(E[i][j - 1]) == get<1>(E[i][j])){
            get<2>(E2.back()) = f(get<2>(E2.back()), get<2>(E[i][j]));
          } else {
            E2.push_back(E[i][j]);
          }
        }
        E[i] = E2;
      }
    }
    P.resize(N + 2);
    P[0] = 0;
    P[1] = 1;
    for (int i = 0; i < N - 1; i++){
      P[i + 2] = P[i + 1] + E[i].size() - 1;
    }
    P[N + 1] = P[N] + 1;
    vector<vector<int>> idL(N), idR(N);
    for (int i = 0; i < N; i++){
      idL[i] = vector<int>(Y[i].size() + 1);
      idR[i] = vector<int>(Y[i].size() + 1);
    }
    idL[0] = vector<int>(Y[0].size() + 1, 0);
    for (int i = 0; i < N - 1; i++){
      int cnt = E[i].size();
      for (int j = 0; j < cnt - 1; j++){
        for (int k = get<0>(E[i][j]) + 1; k <= get<0>(E[i][j + 1]); k++){
          idR[i][k] = P[i + 1] + j;
        }
        for (int k = get<1>(E[i][j]) + 1; k <= get<1>(E[i][j + 1]); k++){
          idL[i + 1][k] = P[i + 1] + j;
        }
      }
    }
    idR[N - 1] = vector<int>(Y[N - 1].size() + 1, P[N]);
    int V = P[N + 1];
    vector<vector<pair<T, int>>> G(V);
    for (int i = 0; i <= N; i++){
      for (int j = P[i]; j < P[i + 1] - 1; j++){
        G[j].push_back(make_pair(get<2>(E[i - 1][j - P[i] + 1]), j + 1));
        G[j + 1].push_back(make_pair(get<2>(E[i - 1][j - P[i] + 1]), j));
      }
    }
    for (int i = 0; i < N; i++){
      for (int j = 0; j <= Y[i].size(); j++){
        G[idL[i][j]].push_back(make_pair(e, idR[i][j]));
        G[idR[i][j]].push_back(make_pair(e, idL[i][j]));
      }
    }
    return G;
  }
  int find(point p){
    int x = lower_bound(X.begin(), X.end(), p.x - eps) - X.begin();
    if (x == 0){
      return 0;
    } else if (x == N){
      return P[N];
    } else {
      int ans = P[x];
      for (int i = 1; i < E[x - 1].size() - 1; i++){
        point A(X[x - 1], Y[x - 1][get<0>(E[x - 1][i])]);
        point B(X[x], Y[x][get<1>(E[x - 1][i])]);
        if (cross(B - A, p - A) > 0){
          ans++;
        }
      }
      return ans;
    }
  }
  vector<double> area(){
    int V = P[N + 1];
    vector<double> ans(V, -1);
    for (int i = 0; i < N - 1; i++){
      for(int j = 1; j < E[i].size() - 2; j++){
        double l = Y[i][get<0>(E[i][j + 1])] - Y[i][get<0>(E[i][j])];
        double r = Y[i + 1][get<1>(E[i][j + 1])] - Y[i + 1][get<1>(E[i][j])];
        ans[P[i + 1] + j] = (l + r) * (X[i + 1] - X[i]) / 2;
      }
    }
    return ans;
  }
  vector<double> length(){
    int V = P[N + 1];
    vector<double> ans(V, -1);
    for (int i = 0; i < N - 1; i++){
      for(int j = 1; j < E[i].size() - 1; j++){
        double l = Y[i][get<0>(E[i][j])];
        double r = Y[i + 1][get<1>(E[i][j])];
        ans[P[i + 1] + j - 1] = hypot(X[i + 1] - X[i], r - l);
      }
    }
    return ans;
  }
};
int main(){
  cout << fixed << setprecision(20);
  while (true){
    int N, R;
    cin >> N >> R;
    if (N == 0 && R == 0){
      break;
    }
    vector<point> P(N + 1);
    for (int i = 0; i < N; i++){
      cin >> P[i].x >> P[i].y;
    }
    P[N] = P[0];
    vector<vector<point>> P2(N, vector<point>(6));
    for (int i = 0; i < N; i++){
      point A = P[i], B = P[i + 1];
      if (A.x + A.y > B.x + B.y){
        swap(A, B);
      }
      if (A.x - A.y > B.x - B.y){
        P2[i][0] = point(A.x - R, A.y);
        P2[i][1] = point(B.x - R, B.y);
        P2[i][2] = point(B.x, B.y + R);
        P2[i][3] = point(B.x + R, B.y);
        P2[i][4] = point(A.x + R, A.y);
        P2[i][5] = point(A.x, A.y - R);
      } else {
        P2[i][0] = point(A.x - R, A.y);
        P2[i][1] = point(A.x, A.y + R);
        P2[i][2] = point(B.x, B.y + R);
        P2[i][3] = point(B.x + R, B.y);
        P2[i][4] = point(B.x, B.y - R);
        P2[i][5] = point(A.x, A.y - R);
      }
    }
    for (int i = 0; i <= N; i++){
      P[i] = rotate(P[i], 1);
    }
    for (int i = 0; i < N; i++){
      for (int j = 0; j < 6; j++){
        P2[i][j] = rotate(P2[i][j], 1);
      }
    }
    slab_decomposition<__int128_t> D(bit_xor<__int128_t>(), 0);
    for (int i = 0; i < N; i++){
      D.add_segment(line(P[i], P[i + 1]), 1);
    }
    for (int i = 0; i < N; i++){
      for (int j = 0; j < 6; j++){
        D.add_segment(line(P2[i][j], P2[i][(j + 1) % 6]), (__int128_t) 1 << (i + 1));
      }
    }
    vector<vector<pair<__int128_t, int>>> E = D.build();
    int V = E.size();
    vector<__int128_t> A(V, -1);
    A[0] = 0;
    queue<int> Q;
    Q.push(0);
    while (!Q.empty()){
      int v = Q.front();
      Q.pop();
      for (pair<__int128_t, int> e : E[v]){
        int w = e.second;
        if (A[w] == -1){
          A[w] = A[v] ^ e.first;
          Q.push(w);
        }
      }
    }
    vector<double> L = D.length();
    double ans = 0;
    for (int i = 0; i < V - 1; i++){
      if (A[i] != 1 && A[i + 1] == 1 || A[i] == 1 && A[i + 1] != 1){
        ans += L[i];
      }
    }
    cout << ans << endl;
  }
}