互いに交差しない多角形のリストが与えられる.適当な点が与えられたとき,それがどの多角形に含まれるかを判定する.
この問題の本来の形は,DCEL で与えられた平面地図上の点位置決定であるが,ここでは入力しやすい形ということで多角形のリストを与えることにしている.教科書的アプローチは台形地図であるが,これを実装するにはそれなりに複雑なアルゴリズムが必要となってしまう.
そこでここでは "Computational Geometry, Algorithms and Applications" で「これは実用的でない」といわれているスラブ地図に基づくアルゴリズムを実装した.なぜ実用的でないかというと,前処理の結果空間計算量が O(n^2) となるからである.
領域の頂点の数を n として,
境界上が 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> ®ions) : 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; } };
Last Modified: 2006.12.11 08:46:08.