点位置決定 (スラブ地図)

説明 (言い訳)

互いに交差しない多角形のリストが与えられる．適当な点が与えられたとき，それがどの多角形に含まれるかを判定する．

この問題の本来の形は，DCEL で与えられた平面地図上の点位置決定であるが，ここでは入力しやすい形ということで多角形のリストを与えることにしている．教科書的アプローチは台形地図であるが，これを実装するにはそれなりに複雑なアルゴリズムが必要となってしまう．

そこでここでは "Computational Geometry, Algorithms and Applications" で「これは実用的でない」といわれているスラブ地図に基づくアルゴリズムを実装した．なぜ実用的でないかというと，前処理の結果空間計算量が O(n^2) となるからである．

計算量

領域の頂点の数を n として，

• 前処理: 時間，空間ともに O(n^2)．
• クエリー: O(log n)．

使い方

point_location(const vector<polygon> &regions)

互いに交差しない多角形のリスト regions を与え，検索しやすい形に前処理する．

int locate(const point &p)

p を含む多角形のインデックスを返す．含む多角形が無い場合は -1 を返す．p が境界上にある場合は don't care．

境界上が don't care となる実装をしているので，あらかじめ多角形を摂動変形しておくことを推奨．このとき変形後の多角形がオーバーラップしないことは問題設定から考慮すること．

TODO

• どこかで Verify する．
• 境界が don't care にならない実装を考える．これは境界共有する多角形への対処にも繋がる．
• 行数の短縮．
• 普通の台形地図を実装したいなあ．
• 実行速度の加速．

ソースコード

#define curr(P,i) P[i]
#define next(P,i) P[(i+1)%P.size()]
#define prev(P,i) P[(i+P.size()-1)%P.size()]
struct point_location {
  vector<polygon> regions;
  point_location(const vector<polygon> &regions) :
    regions(regions) { compile(); }
  struct cut {
    double x, y, a;
    int d;
    cut(double x, double y, double a, int d) :
      x(x), y(y), a(a), d(d) { }
    bool operator < (const cut &c) const {
      return y != c.y ? y < c.y : a != c.a ? a < c.a : d > c.d;
    }
    bool operator == (const cut &c) const {
      return y == c.y && a == c.a;
    }
  };
  vector<double> xs;
  vector< vector<cut> > ys;

  void compile() {
    for (int i = 0; i < regions.size(); ++i)
      for (int j = 0; j < regions[i].size(); ++j)
        xs.push_back(real(regions[i][j]));
    sort(xs.begin(), xs.end());
    xs.erase(unique(xs.begin(), xs.end()), xs.end());
    ys.resize(xs.size());

    for (int k = 0; k < xs.size(); ++k) {
      for (int i = 0; i < regions.size(); ++i) {
        for (int j = 0; j < regions[i].size(); ++j) {
          segment seg( curr(regions[i],j), next(regions[i],j) );
          int dominant = i;
          if (real(seg[0]) > real(seg[1])) {
            swap(seg[0], seg[1]);
            dominant = -1;
          }
          if (xs[k] < real(seg[0]) || real(seg[1]) <= xs[k]) continue;

          double a = imag(seg[1]-seg[0])/real(seg[1]-seg[0]);
          double y = imag(seg[0]) + a * (xs[k]-real(seg[0]));
          ys[k].push_back( cut(xs[k], y, a, dominant) );
        }
      }
      sort(ys[k].begin(), ys[k].end());
      ys[k].erase(unique(ys[k].begin(), ys[k].end()), ys[k].end());
    }
  }
  struct at_x {
    double X;
    at_x(double X) : X(X) { }
    bool operator () (double y, const cut &c) const {
      return y < c.y + c.a * (X - c.x);
    }
  };
  int locate(const point &p) {
    double x = real(p), y = imag(p);
    int i = distance(xs.begin(),
        upper_bound(xs.begin(), xs.end(), x)) - 1;
    if (i < 0 || i >= xs.size()) return -1;
    int j = distance(ys[i].begin(),
        upper_bound(ys[i].begin(), ys[i].end(), y, at_x(x))) - 1;
    if (j < 0 || j >= ys[i].size()) return -1;
    return ys[i][j].d;
  }
};

Verified