2D ポイントがポリゴン内にあるかどうかを確認するにはどうすればよいですか?
https://stackoverflow.com/questions/217578
-
03-07-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
Russian質問
を作成しようとしています 速い ポリゴン内の 2D ポイント アルゴリズム。ヒット テストで使用します (例: Polygon.contains(p:Point)）。効果的なテクニックの提案をいただければ幸いです。
解決
グラフィックスの場合、私は整数を好みません。多くのシステムは UI ペイントに整数を使用します (結局、ピクセルは整数です) が、たとえば macOS ではすべてに float が使用されます。macOS はポイントのみを認識し、ポイントは 1 つのピクセルに変換できますが、モニターの解像度によっては、他のものに変換される場合があります。Retina スクリーンでは、半分のポイント (0.5/0.5) がピクセルです。それでも、macOS UI が他の UI に比べて大幅に遅いことに私は気づきませんでした。結局のところ、3D API (OpenGL または Direct3D) は float でも動作し、最新のグラフィックス ライブラリは GPU アクセラレーションを利用することが非常に多いのです。
さて、スピードが最大の関心事だと言いましたが、わかりました、スピードを追求しましょう。高度なアルゴリズムを実行する前に、まず簡単なテストを実行します。を作成します 軸に合わせたバウンディングボックス ポリゴンの周りに。これは非常に簡単かつ高速で、すでに多くの計算を安全に行うことができます。それはどのように機能するのでしょうか?多角形のすべての点を反復処理し、X と Y の最小値/最大値を見つけます。
例えば。あなたはポイントを持っています (9/1), (4/3), (2/7), (8/2), (3/6). 。これは、Xmin が 2、Xmax が 9、Ymin が 1、Ymax が 7 であることを意味します。2 つのエッジ (2/1) と (9/7) を持つ長方形の外側の点を多角形の中に入れることはできません。
// p is your point, p.x is the x coord, p.y is the y coord
if (p.x < Xmin || p.x > Xmax || p.y < Ymin || p.y > Ymax) {
// Definitely not within the polygon!
}
これは、任意のポイントに対して実行される最初のテストです。ご覧のとおり、このテストは非常に高速ですが、非常に粗いテストでもあります。外接する四角形内の点を処理するには、より洗練されたアルゴリズムが必要です。これを計算するにはいくつかの方法があります。どちらの方法が機能するかは、ポリゴンに穴がある可能性があるか、それとも常にソリッドであるかによって決まります。ソリッドの例 (凸面 1 つ、凹面 1 つ) を次に示します。
そして、これが穴のあるものです。
緑色のものは真ん中に穴が開いています！
上記の 3 つのケースすべてを処理でき、それでもかなり高速な、最も簡単なアルゴリズムは次のとおりです。 レイキャスティング. 。アルゴリズムの考え方は非常にシンプルです。ポリゴンの外側の任意の場所からポイントまで仮想レイを描画し、それがポリゴンの側面に当たる頻度を数えます。ヒット数が偶数の場合はポリゴンの外側にあり、奇数の場合はポリゴンの内側にあります。
の 巻数アルゴリズム 代替手段として、多角形の線に非常に近い点の場合はより正確になりますが、速度も大幅に遅くなります。浮動小数点精度の制限と丸めの問題のため、ポリゴンの辺に近すぎる点ではレイ キャスティングが失敗することがありますが、実際にはそれはほとんど問題ではありません。点が辺の近くにある場合、多くの場合視覚的に不可能であるためです。ビューアは、すでに屋内にいるのか、まだ屋外にいるのかを認識します。
上記の境界ボックスがまだ残っているのを覚えていますか?境界ボックスの外側の点を選択し、それを光線の開始点として使用するだけです。例えば。ポイント (Xmin - e/p.y) 確実にポリゴンの外側にあります。
しかし、何ですか e?良い、 e (実際にはイプシロン) 境界ボックスにいくらかの値を与えます パディング. 。先ほども述べたように、ポリゴン ラインに近づきすぎるとレイ トレーシングは失敗します。境界ボックスは多角形と等しい可能性があるため (多角形が軸に揃えられた長方形の場合、境界ボックスは多角形自体と等しい!)、これを安全にするためにある程度のパディングが必要です。それだけです。どれくらいの大きさを選ぶべきか e?大きすぎません。それは、描画に使用する座標系スケールによって異なります。ピクセル ステップ幅が 1.0 の場合は、1.0 を選択します (0.1 でも同様に機能します)。
開始座標と終了座標を含む光線が得られたので、問題は「」から変わります。多角形内の点です" に "光線はどのくらいの頻度でポリゴンの辺と交差しますか」。したがって、以前のように多角形の点だけを操作することはできず、実際の辺が必要になります。辺は常に 2 点によって定義されます。
side 1: (X1/Y1)-(X2/Y2)
side 2: (X2/Y2)-(X3/Y3)
side 3: (X3/Y3)-(X4/Y4)
:
すべての面に対して光線をテストする必要があります。光線をベクトル、すべての辺をベクトルと考えてください。光線は各面に 1 回ずつ当たるか、まったく当たらない必要があります。同じ面に二度攻撃することはできません。2D 空間内の 2 本の線は、平行でない限り、常に 1 回だけ交差します。平行な場合、交差することはありません。ただし、ベクトルの長さは限られているため、2 つのベクトルは平行にならず、短すぎて互いに交わらない可能性があります。
// Test the ray against all sides
int intersections = 0;
for (side = 0; side < numberOfSides; side++) {
// Test if current side intersects with ray.
// If yes, intersections++;
}
if ((intersections & 1) == 1) {
// Inside of polygon
} else {
// Outside of polygon
}
ここまではうまくいきましたが、2 つのベクトルが交差するかどうかをどのようにテストするのでしょうか?以下にいくつかの C コード (テストされていません) を示します。これでうまくいくはずです。
#define NO 0
#define YES 1
#define COLLINEAR 2
int areIntersecting(
float v1x1, float v1y1, float v1x2, float v1y2,
float v2x1, float v2y1, float v2x2, float v2y2
) {
float d1, d2;
float a1, a2, b1, b2, c1, c2;
// Convert vector 1 to a line (line 1) of infinite length.
// We want the line in linear equation standard form: A*x + B*y + C = 0
// See: http://en.wikipedia.org/wiki/Linear_equation
a1 = v1y2 - v1y1;
b1 = v1x1 - v1x2;
c1 = (v1x2 * v1y1) - (v1x1 * v1y2);
// Every point (x,y), that solves the equation above, is on the line,
// every point that does not solve it, is not. The equation will have a
// positive result if it is on one side of the line and a negative one
// if is on the other side of it. We insert (x1,y1) and (x2,y2) of vector
// 2 into the equation above.
d1 = (a1 * v2x1) + (b1 * v2y1) + c1;
d2 = (a1 * v2x2) + (b1 * v2y2) + c1;
// If d1 and d2 both have the same sign, they are both on the same side
// of our line 1 and in that case no intersection is possible. Careful,
// 0 is a special case, that's why we don't test ">=" and "<=",
// but "<" and ">".
if (d1 > 0 && d2 > 0) return NO;
if (d1 < 0 && d2 < 0) return NO;
// The fact that vector 2 intersected the infinite line 1 above doesn't
// mean it also intersects the vector 1. Vector 1 is only a subset of that
// infinite line 1, so it may have intersected that line before the vector
// started or after it ended. To know for sure, we have to repeat the
// the same test the other way round. We start by calculating the
// infinite line 2 in linear equation standard form.
a2 = v2y2 - v2y1;
b2 = v2x1 - v2x2;
c2 = (v2x2 * v2y1) - (v2x1 * v2y2);
// Calculate d1 and d2 again, this time using points of vector 1.
d1 = (a2 * v1x1) + (b2 * v1y1) + c2;
d2 = (a2 * v1x2) + (b2 * v1y2) + c2;
// Again, if both have the same sign (and neither one is 0),
// no intersection is possible.
if (d1 > 0 && d2 > 0) return NO;
if (d1 < 0 && d2 < 0) return NO;
// If we get here, only two possibilities are left. Either the two
// vectors intersect in exactly one point or they are collinear, which
// means they intersect in any number of points from zero to infinite.
if ((a1 * b2) - (a2 * b1) == 0.0f) return COLLINEAR;
// If they are not collinear, they must intersect in exactly one point.
return YES;
}
入力値は次のとおりです。 2つのエンドポイント ベクトル 1 (v1x1/v1y1 そして v1x2/v1y2) とベクトル 2 (v2x1/v2y1 そして v2x2/v2y2）。つまり、2 つのベクトル、4 つの点、8 つの座標があります。 YES そして NO は明らかです。 YES 交差点が増えて、 NO 何もしません。
コリニアはどうでしょうか？これは、両方のベクトルが同じ無限の線上にあり、位置と長さに応じて、まったく交差しないか、無限の数の点で交差することを意味します。このケースをどのように処理するか完全にはわかりませんが、どちらにしても交差点としてカウントしません。まあ、浮動小数点の丸め誤差があるため、このケースは実際にはかなりまれです。より良いコードはおそらくテストしないでしょう == 0.0f しかし、代わりに次のようなもののために < epsilon, ここで、イプシロンはかなり小さな数です。
より多くのポイントをテストする必要がある場合は、多角形の辺の一次方程式の標準形式をメモリに保存しておくと、毎回再計算する必要がなく、確実に全体を少し高速化できます。これにより、ポリゴンの辺ごとに 3 つの浮動小数点値をメモリに保存する代わりに、テストごとに 2 つの浮動小数点乗算と 3 つの浮動小数点減算が節約されます。これは典型的なメモリと計算時間のトレードオフです。
最後になりましたが、重要なことです：3D ハードウェアを使用して問題を解決できる場合は、興味深い代替手段があります。すべての作業を GPU に任せてください。画面外にペイント サーフェスを作成します。黒で完全に塗りつぶします。次に、OpenGL または Direct3D でポリゴン (または、ポイントがいずれかのポリゴン内にあるかどうかをテストしたい場合はすべてのポリゴンでも) をペイントし、そのポリゴンを別の色で塗りつぶします。色、例:白。点がポリゴン内にあるかどうかを確認するには、描画面からこの点の色を取得します。これは単なる O(1) メモリ フェッチです。
もちろん、この方法は描画面を大きくする必要がない場合にのみ使用できます。GPU メモリに収まらない場合、この方法は CPU で実行するよりも遅くなります。それが巨大である必要があり、GPU が最新のシェーダーをサポートしている場合でも、上記のレイ キャスティングを GPU シェーダーとして実装することで GPU を使用できます。これは絶対に可能です。より多くのポリゴンやテストするポイントの場合、これは効果的です。一部の GPU では 64 ～ 256 のポイントを並行してテストできることを考慮してください。ただし、CPU から GPU へのデータ転送、およびその逆のデータ転送は常にコストがかかることに注意してください。そのため、ポイントまたはポリゴンのいずれかが動的で頻繁に変化する、いくつかの単純なポリゴンに対していくつかのポイントをテストするだけの場合、GPU アプローチが利益をもたらすことはほとんどありません。オフ。
他のヒント
次のコードが最良の解決策だと思います（こちら）：
int pnpoly(int nvert, float *vertx, float *verty, float testx, float testy)
{
int i, j, c = 0;
for (i = 0, j = nvert-1; i < nvert; j = i++) {
if ( ((verty[i]>testy) != (verty[j]>testy)) &&
(testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) )
c = !c;
}
return c;
}
引数
- nvert ：ポリゴン内の頂点の数。最後に最初の頂点を繰り返すかどうかは、上記の記事で説明されています。
- vertx、verty ：ポリゴンの頂点のx座標とy座標を含む配列。
- testx、testy ：テストポイントのX座標とY座標。
短くて効率的で、凸面と凹面の両方のポリゴンで機能します。前に提案したように、最初に境界矩形をチェックし、多角形の穴を個別に処理する必要があります。
この背後にある考え方は非常に単純です。著者は次のように説明しています：
テストポイントから半無限光線を水平方向に（xを増やし、yを固定して）実行し、それが交差するエッジの数をカウントします。交差するたびに、光線は内側と外側を切り替えます。これは、ヨルダン曲線定理と呼ばれます。
変数cは、水平光線がエッジを横切るたびに0から1および1から0に切り替わります。したがって、基本的には、交差するエッジの数が偶数か奇数かを追跡しています。 0は偶数、1は奇数を意味します。
こちらは、 nirgによる回答のC＃バージョンです。このRPI教授。そのRPIソースのコードを使用するには属性が必要です。
境界ボックスチェックが上部に追加されました。ただし、James Brownが指摘しているように、メインコードはバウンディングボックスチェック自体とほぼ同じ速度であるため、チェックするポイントのほとんどがバウンディングボックス内にある場合、バウンディングボックスチェックは実際に全体の動作を遅くする可能性があります。そのため、バウンディングボックスをチェックアウトしたままにすることもできますが、あまり頻繁に形状が変化しない場合は、ポリゴンのバウンディングボックスを事前に計算することもできます。
public bool IsPointInPolygon( Point p, Point[] polygon )
{
double minX = polygon[ 0 ].X;
double maxX = polygon[ 0 ].X;
double minY = polygon[ 0 ].Y;
double maxY = polygon[ 0 ].Y;
for ( int i = 1 ; i < polygon.Length ; i++ )
{
Point q = polygon[ i ];
minX = Math.Min( q.X, minX );
maxX = Math.Max( q.X, maxX );
minY = Math.Min( q.Y, minY );
maxY = Math.Max( q.Y, maxY );
}
if ( p.X < minX || p.X > maxX || p.Y < minY || p.Y > maxY )
{
return false;
}
// http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
bool inside = false;
for ( int i = 0, j = polygon.Length - 1 ; i < polygon.Length ; j = i++ )
{
if ( ( polygon[ i ].Y > p.Y ) != ( polygon[ j ].Y > p.Y ) &&
p.X < ( polygon[ j ].X - polygon[ i ].X ) * ( p.Y - polygon[ i ].Y ) / ( polygon[ j ].Y - polygon[ i ].Y ) + polygon[ i ].X )
{
inside = !inside;
}
}
return inside;
}
これは、Nirgのアプローチに基づいたM. Katzによる回答のJavaScriptバリアントです。
function pointIsInPoly(p, polygon) {
var isInside = false;
var minX = polygon[0].x, maxX = polygon[0].x;
var minY = polygon[0].y, maxY = polygon[0].y;
for (var n = 1; n < polygon.length; n++) {
var q = polygon[n];
minX = Math.min(q.x, minX);
maxX = Math.max(q.x, maxX);
minY = Math.min(q.y, minY);
maxY = Math.max(q.y, maxY);
}
if (p.x < minX || p.x > maxX || p.y < minY || p.y > maxY) {
return false;
}
var i = 0, j = polygon.length - 1;
for (i, j; i < polygon.length; j = i++) {
if ( (polygon[i].y > p.y) != (polygon[j].y > p.y) &&
p.x < (polygon[j].x - polygon[i].x) * (p.y - polygon[i].y) / (polygon[j].y - polygon[i].y) + polygon[i].x ) {
isInside = !isInside;
}
}
return isInside;
}
点pと各多角形の頂点の間の角度の方向の合計を計算します。指向角の合計が360度の場合、ポイントは内側にあります。合計が0の場合、ポイントは外側です。
この方法はより堅牢で、数値の精度に依存しにくいため、この方法のほうが好きです。
交差点の数の計算中に頂点を「ヒット」できるため、交差点の数の均一性を計算する方法は制限されます。
編集：ところで、この方法は凹面と凸面のポリゴンで機能します。
編集：最近、トピックに関するウィキペディアの記事を見つけました。
boboboboが引用した Eric Hainesの記事は非常に優れています。特に興味深いのは、アルゴリズムのパフォーマンスを比較した表です。角度の加算方法は他の方法と比べて本当に悪いです。また興味深いのは、ルックアップグリッドを使用してポリゴンを<！> quot; in <！> quotにさらに細分化するような最適化です。および<！> quot; out <！> quot;セクターは、<！> gtのあるポリゴンでも非常に高速にテストを実行できます。 1000辺。
とにかく、まだ早いですが、私の投票は<！> quot; crossings <！> quot;に行きます。これは、Meckiが説明しているとおりです。しかし、私はそれを最も簡潔に見つけました David Bourkeによって記述および成文化されました。実際の三角法は不要であり、凸面と凹面で機能し、辺の数が増えても十分に機能します。
ところで、エリックヘインズの興味のある記事からのパフォーマンステーブルの1つは、ランダムポリゴンでのテストです。
number of edges per polygon
3 4 10 100 1000
MacMartin 2.9 3.2 5.9 50.6 485
Crossings 3.1 3.4 6.8 60.0 624
Triangle Fan+edge sort 1.1 1.8 6.5 77.6 787
Triangle Fan 1.2 2.1 7.3 85.4 865
Barycentric 2.1 3.8 13.8 160.7 1665
Angle Summation 56.2 70.4 153.6 1403.8 14693
Grid (100x100) 1.5 1.5 1.6 2.1 9.8
Grid (20x20) 1.7 1.7 1.9 5.7 42.2
Bins (100) 1.8 1.9 2.7 15.1 117
Bins (20) 2.1 2.2 3.7 26.3 278
この質問はとても興味深いです。この投稿の他の回答とは異なる別の実行可能なアイデアがあります。
アイデアは、角度の合計を使用して、ターゲットが内側か外側かを判断することです。ターゲットがエリア内にある場合、ターゲットと2つの境界点ごとの角度の合計は360になります。ターゲットが外部にある場合、合計は360にはなりません。角度には方向があります。角度が後方に行く場合、角度は負の角度です。これは、ワインディング番号を計算するようなものです。
このイメージを参照して、アイデアの基本を理解してください。
私のアルゴリズムでは、時計回りが正の方向であると想定しています。潜在的な入力は次のとおりです。
[[-122.402015, 48.225216], [-117.032049, 48.999931], [-116.919132, 45.995175], [-124.079107, 46.267259], [-124.717175, 48.377557], [-122.92315, 47.047963], [-122.402015, 48.225216]]
以下は、アイデアを実装するpythonコードです。
def isInside(self, border, target):
degree = 0
for i in range(len(border) - 1):
a = border[i]
b = border[i + 1]
# calculate distance of vector
A = getDistance(a[0], a[1], b[0], b[1]);
B = getDistance(target[0], target[1], a[0], a[1])
C = getDistance(target[0], target[1], b[0], b[1])
# calculate direction of vector
ta_x = a[0] - target[0]
ta_y = a[1] - target[1]
tb_x = b[0] - target[0]
tb_y = b[1] - target[1]
cross = tb_y * ta_x - tb_x * ta_y
clockwise = cross < 0
# calculate sum of angles
if(clockwise):
degree = degree + math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C)))
else:
degree = degree - math.degrees(math.acos((B * B + C * C - A * A) / (2.0 * B * C)))
if(abs(round(degree) - 360) <= 3):
return True
return False
answer by nirg の迅速なバージョン：
extension CGPoint {
func isInsidePolygon(vertices: [CGPoint]) -> Bool {
guard !vertices.isEmpty else { return false }
var j = vertices.last!, c = false
for i in vertices {
let a = (i.y > y) != (j.y > y)
let b = (x < (j.x - i.x) * (y - i.y) / (j.y - i.y) + i.x)
if a && b { c = !c }
j = i
}
return c
}
}
Michael Stonebraker の下で研究者だったときに、この作業をやり直しました。 Ingres 、 PostgreSQL など
超高速なので、最初にバウンディングボックスを作成することが最速の方法であることがわかりました。境界ボックスの外側にある場合は、外側にあります。それ以外の場合は、より難しい作業を行います...
優れたアルゴリズムが必要な場合は、ジオワーク用のオープンソースプロジェクトPostgreSQLソースコードをご覧ください...
指摘したいのですが、右利き対左利きについての洞察は得られませんでした（<！> quot; inside <！> quot; vs <！> quot; outside <！> quot;問題としても表現できます。 。
更新
BKBのリンクは、かなりの数の合理的なアルゴリズムを提供しました。私は地球科学の問題に取り組んでいたので、緯度/経度で動作する解決策が必要でしたが、それには利き手の特有の問題があります-小さい領域の内側の領域ですか、大きい領域ですか？答えは、<！> quot; direction <！> quot;頂点の問題-それは左利きまたは右利きのいずれかであり、この方法では、どちらかのエリアを<！> quot; inside <！> quotとして示すことができます。任意のポリゴン。そのため、私の作業では、そのページに列挙されているソリューション3を使用しました。
さらに、私の仕事では、<！> quot;行の<！> quot;に個別の関数を使用しました。テスト。
...誰かが尋ねたので：頂点の数が特定の数を超えたときにバウンディングボックステストが最適であることがわかりました-必要に応じて、より長いテストを行う前に非常に迅速なテストを行います...単に最大のx、最小のx、最大のy、最小のyを取り、それらを組み合わせてボックスの4つのポイントを作成します...
次のヒント：もう少し洗練された<！> quot; light-dimming <！> quot;平面上のすべての正の点でグリッド空間で計算し、次に<！> quot; real <！> quotに再投影します。経度/緯度。これにより、経度180の線を横切ったときや極地域を処理したときに、ラップする可能性のあるエラーを回避できます。うまくいきました！
Nirgが投稿し、boboboboが編集したソリューションに本当に似ています。 JavaScriptを使いやすくし、もう少し読みやすくしました：
function insidePoly(poly, pointx, pointy) {
var i, j;
var inside = false;
for (i = 0, j = poly.length - 1; i < poly.length; j = i++) {
if(((poly[i].y > pointy) != (poly[j].y > pointy)) && (pointx < (poly[j].x-poly[i].x) * (pointy-poly[i].y) / (poly[j].y-poly[i].y) + poly[i].x) ) inside = !inside;
}
return inside;
}
David Segondの答えはほぼ標準的な一般的な答えであり、Richard Tの答えは最も一般的な最適化ですが、他にもいくつかあります。他の強力な最適化は、一般的ではないソリューションに基づいています。たとえば、同じポリゴンを多数のポイントでチェックする場合、非常に高速なTIN検索アルゴリズムが多数あるため、ポリゴンを三角形分割すると速度が大幅に向上します。もう1つは、画面表示などの低解像度でポリゴンとポイントが制限された平面上にある場合、ポリゴンを特定の色でメモリマップ表示バッファーにペイントし、特定のピクセルの色をチェックして、それが存在するかどうかを確認することですポリゴン内。
多くの最適化と同様に、これらは一般的なケースではなく特定のケースに基づいており、1回の使用ではなく償却時間に基づいたメリットをもたらします。
この分野で働いていると、Joeseph O'Rourkesの「Computation Geometry in C」ISBN 0-521-44034-3が非常に役立つことがわかりました。
簡単な解決策は、こちら
ポリゴンが CONVEX の場合、より良いアプローチがあるかもしれません。無限の線の集まりとして多角形を見てください。スペースを2つに分割する各行。すべてのポイントについて、それが線の片側にあるのか、それとも反対側にあるのかは簡単です。ポイントがすべてのラインの同じ側にある場合、それはポリゴン内にあります。
これは古いことに気づきましたが、興味がある人のために、ココアに実装されたレイキャスティングアルゴリズムを以下に示します。それが物事を行うための最も効率的な方法であるかどうかはわかりませんが、誰かを助けるかもしれません。
- (BOOL)shape:(NSBezierPath *)path containsPoint:(NSPoint)point
{
NSBezierPath *currentPath = [path bezierPathByFlatteningPath];
BOOL result;
float aggregateX = 0; //I use these to calculate the centroid of the shape
float aggregateY = 0;
NSPoint firstPoint[1];
[currentPath elementAtIndex:0 associatedPoints:firstPoint];
float olderX = firstPoint[0].x;
float olderY = firstPoint[0].y;
NSPoint interPoint;
int noOfIntersections = 0;
for (int n = 0; n < [currentPath elementCount]; n++) {
NSPoint points[1];
[currentPath elementAtIndex:n associatedPoints:points];
aggregateX += points[0].x;
aggregateY += points[0].y;
}
for (int n = 0; n < [currentPath elementCount]; n++) {
NSPoint points[1];
[currentPath elementAtIndex:n associatedPoints:points];
//line equations in Ax + By = C form
float _A_FOO = (aggregateY/[currentPath elementCount]) - point.y;
float _B_FOO = point.x - (aggregateX/[currentPath elementCount]);
float _C_FOO = (_A_FOO * point.x) + (_B_FOO * point.y);
float _A_BAR = olderY - points[0].y;
float _B_BAR = points[0].x - olderX;
float _C_BAR = (_A_BAR * olderX) + (_B_BAR * olderY);
float det = (_A_FOO * _B_BAR) - (_A_BAR * _B_FOO);
if (det != 0) {
//intersection points with the edges
float xIntersectionPoint = ((_B_BAR * _C_FOO) - (_B_FOO * _C_BAR)) / det;
float yIntersectionPoint = ((_A_FOO * _C_BAR) - (_A_BAR * _C_FOO)) / det;
interPoint = NSMakePoint(xIntersectionPoint, yIntersectionPoint);
if (olderX <= points[0].x) {
//doesn't matter in which direction the ray goes, so I send it right-ward.
if ((interPoint.x >= olderX && interPoint.x <= points[0].x) && (interPoint.x > point.x)) {
noOfIntersections++;
}
} else {
if ((interPoint.x >= points[0].x && interPoint.x <= olderX) && (interPoint.x > point.x)) {
noOfIntersections++;
}
}
}
olderX = points[0].x;
olderY = points[0].y;
}
if (noOfIntersections % 2 == 0) {
result = FALSE;
} else {
result = TRUE;
}
return result;
}
ポイントをテストするためのサンプルメソッドを使用したnirgの答えのObj-Cバージョン。 Nirgの答えは私にとってはうまくいきました。
- (BOOL)isPointInPolygon:(NSArray *)vertices point:(CGPoint)test {
NSUInteger nvert = [vertices count];
NSInteger i, j, c = 0;
CGPoint verti, vertj;
for (i = 0, j = nvert-1; i < nvert; j = i++) {
verti = [(NSValue *)[vertices objectAtIndex:i] CGPointValue];
vertj = [(NSValue *)[vertices objectAtIndex:j] CGPointValue];
if (( (verti.y > test.y) != (vertj.y > test.y) ) &&
( test.x < ( vertj.x - verti.x ) * ( test.y - verti.y ) / ( vertj.y - verti.y ) + verti.x) )
c = !c;
}
return (c ? YES : NO);
}
- (void)testPoint {
NSArray *polygonVertices = [NSArray arrayWithObjects:
[NSValue valueWithCGPoint:CGPointMake(13.5, 41.5)],
[NSValue valueWithCGPoint:CGPointMake(42.5, 56.5)],
[NSValue valueWithCGPoint:CGPointMake(39.5, 69.5)],
[NSValue valueWithCGPoint:CGPointMake(42.5, 84.5)],
[NSValue valueWithCGPoint:CGPointMake(13.5, 100.0)],
[NSValue valueWithCGPoint:CGPointMake(6.0, 70.5)],
nil
];
CGPoint tappedPoint = CGPointMake(23.0, 70.0);
if ([self isPointInPolygon:polygonVertices point:tappedPoint]) {
NSLog(@"YES");
} else {
NSLog(@"NO");
}
}
nirgの答えのC＃バージョンはこちらです：コードを共有するだけです。時間を節約できる可能性があります。
public static bool IsPointInPolygon(IList<Point> polygon, Point testPoint) {
bool result = false;
int j = polygon.Count() - 1;
for (int i = 0; i < polygon.Count(); i++) {
if (polygon[i].Y < testPoint.Y && polygon[j].Y >= testPoint.Y || polygon[j].Y < testPoint.Y && polygon[i].Y >= testPoint.Y) {
if (polygon[i].X + (testPoint.Y - polygon[i].Y) / (polygon[j].Y - polygon[i].Y) * (polygon[j].X - polygon[i].X) < testPoint.X) {
result = !result;
}
}
j = i;
}
return result;
}
問題の帰納的定義ほど美しいものはありません。ここでは完全を期すために、プロローグにバージョンがあります。これにより、レイキャスティングの背後にある考えも明確になります。
http：//の単純化アルゴリズムのシミュレーションに基づくwww.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
ヘルパー述語：
exor(A,B):- \+A,B;A,\+B.
in_range(Coordinate,CA,CB) :- exor((CA>Coordinate),(CB>Coordinate)).
inside(false).
inside(_,[_|[]]).
inside(X:Y, [X1:Y1,X2:Y2|R]) :- in_range(Y,Y1,Y2), X > ( ((X2-X1)*(Y-Y1))/(Y2-Y1) + X1),toggle_ray, inside(X:Y, [X2:Y2|R]); inside(X:Y, [X2:Y2|R]).
get_line(_,_,[]).
get_line([XA:YA,XB:YB],[X1:Y1,X2:Y2|R]):- [XA:YA,XB:YB]=[X1:Y1,X2:Y2]; get_line([XA:YA,XB:YB],[X2:Y2|R]).
AとBの2点が与えられた線の方程式（Line（A、B））は次のとおりです。
(YB-YA)
Y - YA = ------- * (X - XA)
(XB-YB)
線の回転方向が 境界は時計回りに、穴は反時計回りに設定されています。 ポイント（X、Y）、つまりテストされたポイントが左側にあるかどうかを確認します 私たちのラインの半平面（それは好みの問題であり、 右側だけでなく、境界線の方向も変更する必要があります その場合）、これは点から右（または左）に光線を投影することです 線との交差点を確認します。投影することを選択しました 水平方向の光線（これも好みの問題です。 同様の制限付きで垂直に行うこともできます）：
(XB-XA)
X < ------- * (Y - YA) + XA
(YB-YA)
ここで、ポイントが左側（または右側）にあるかどうかを知る必要があります 平面全体ではなく、線分のみであるため、 このセグメントのみに検索を制限しますが、これは簡単です。 セグメント内にあるためには、ライン内の1ポイントだけが高くなる 縦軸のYより。これは強い制限なので、 最初にチェックする必要があるため、これらの行のみを最初に取得します この要件を満たし、その位置を確認します。ヨルダン 多角形に投影される光線は、 偶数行。これで完了です。レイを 正しく、それが線と交差するたびに、その状態を切り替えます。 ただし、実装では、次の長さをチェックします。 指定された制限を満たし、決定するソリューションのバッグ それに内なる。ポリゴン内の各ラインに対してこれを行う必要があります。
is_left_half_plane(_,[],[],_).
is_left_half_plane(X:Y,[XA:YA,XB:YB], [[X1:Y1,X2:Y2]|R], Test) :- [XA:YA, XB:YB] = [X1:Y1, X2:Y2], call(Test, X , (((XB - XA) * (Y - YA)) / (YB - YA) + XA));
is_left_half_plane(X:Y, [XA:YA, XB:YB], R, Test).
in_y_range_at_poly(Y,[XA:YA,XB:YB],Polygon) :- get_line([XA:YA,XB:YB],Polygon), in_range(Y,YA,YB).
all_in_range(Coordinate,Polygon,Lines) :- aggregate(bag(Line), in_y_range_at_poly(Coordinate,Line,Polygon), Lines).
traverses_ray(X:Y, Lines, Count) :- aggregate(bag(Line), is_left_half_plane(X:Y, Line, Lines, <), IntersectingLines), length(IntersectingLines, Count).
% This is the entry point predicate
inside_poly(X:Y,Polygon,Answer) :- all_in_range(Y,Polygon,Lines), traverses_ray(X:Y, Lines, Count), (1 is mod(Count,2)->Answer=inside;Answer=outside).
Javaバージョン：
public class Geocode {
private float latitude;
private float longitude;
public Geocode() {
}
public Geocode(float latitude, float longitude) {
this.latitude = latitude;
this.longitude = longitude;
}
public float getLatitude() {
return latitude;
}
public void setLatitude(float latitude) {
this.latitude = latitude;
}
public float getLongitude() {
return longitude;
}
public void setLongitude(float longitude) {
this.longitude = longitude;
}
}
public class GeoPolygon {
private ArrayList<Geocode> points;
public GeoPolygon() {
this.points = new ArrayList<Geocode>();
}
public GeoPolygon(ArrayList<Geocode> points) {
this.points = points;
}
public GeoPolygon add(Geocode geo) {
points.add(geo);
return this;
}
public boolean inside(Geocode geo) {
int i, j;
boolean c = false;
for (i = 0, j = points.size() - 1; i < points.size(); j = i++) {
if (((points.get(i).getLongitude() > geo.getLongitude()) != (points.get(j).getLongitude() > geo.getLongitude())) &&
(geo.getLatitude() < (points.get(j).getLatitude() - points.get(i).getLatitude()) * (geo.getLongitude() - points.get(i).getLongitude()) / (points.get(j).getLongitude() - points.get(i).getLongitude()) + points.get(i).getLatitude()))
c = !c;
}
return c;
}
}
.Netポート：
static void Main(string[] args)
{
Console.Write("Hola");
List<double> vertx = new List<double>();
List<double> verty = new List<double>();
int i, j, c = 0;
vertx.Add(1);
vertx.Add(2);
vertx.Add(1);
vertx.Add(4);
vertx.Add(4);
vertx.Add(1);
verty.Add(1);
verty.Add(2);
verty.Add(4);
verty.Add(4);
verty.Add(1);
verty.Add(1);
int nvert = 6; //Vértices del poligono
double testx = 2;
double testy = 5;
for (i = 0, j = nvert - 1; i < nvert; j = i++)
{
if (((verty[i] > testy) != (verty[j] > testy)) &&
(testx < (vertx[j] - vertx[i]) * (testy - verty[i]) / (verty[j] - verty[i]) + vertx[i]))
c = 1;
}
}
VBAバージョン：
注：ポリゴンがマップ内の領域である場合、緯度/経度はX / Y（緯度= Y、経度= X）ではなくY / X値であることに注意してください。これは、経度が測定値ではなかったときの方法。
クラスモジュール：CPoint
Private pXValue As Double
Private pYValue As Double
'''''X Value Property'''''
Public Property Get X() As Double
X = pXValue
End Property
Public Property Let X(Value As Double)
pXValue = Value
End Property
'''''Y Value Property'''''
Public Property Get Y() As Double
Y = pYValue
End Property
Public Property Let Y(Value As Double)
pYValue = Value
End Property
モジュール：
Public Function isPointInPolygon(p As CPoint, polygon() As CPoint) As Boolean
Dim i As Integer
Dim j As Integer
Dim q As Object
Dim minX As Double
Dim maxX As Double
Dim minY As Double
Dim maxY As Double
minX = polygon(0).X
maxX = polygon(0).X
minY = polygon(0).Y
maxY = polygon(0).Y
For i = 1 To UBound(polygon)
Set q = polygon(i)
minX = vbMin(q.X, minX)
maxX = vbMax(q.X, maxX)
minY = vbMin(q.Y, minY)
maxY = vbMax(q.Y, maxY)
Next i
If p.X < minX Or p.X > maxX Or p.Y < minY Or p.Y > maxY Then
isPointInPolygon = False
Exit Function
End If
' SOURCE: http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
isPointInPolygon = False
i = 0
j = UBound(polygon)
Do While i < UBound(polygon) + 1
If (polygon(i).Y > p.Y) Then
If (polygon(j).Y < p.Y) Then
If p.X < (polygon(j).X - polygon(i).X) * (p.Y - polygon(i).Y) / (polygon(j).Y - polygon(i).Y) + polygon(i).X Then
isPointInPolygon = True
Exit Function
End If
End If
ElseIf (polygon(i).Y < p.Y) Then
If (polygon(j).Y > p.Y) Then
If p.X < (polygon(j).X - polygon(i).X) * (p.Y - polygon(i).Y) / (polygon(j).Y - polygon(i).Y) + polygon(i).X Then
isPointInPolygon = True
Exit Function
End If
End If
End If
j = i
i = i + 1
Loop
End Function
Function vbMax(n1, n2) As Double
vbMax = IIf(n1 > n2, n1, n2)
End Function
Function vbMin(n1, n2) As Double
vbMin = IIf(n1 > n2, n2, n1)
End Function
Sub TestPointInPolygon()
Dim i As Integer
Dim InPolygon As Boolean
' MARKER Object
Dim p As CPoint
Set p = New CPoint
p.X = <ENTER X VALUE HERE>
p.Y = <ENTER Y VALUE HERE>
' POLYGON OBJECT
Dim polygon() As CPoint
ReDim polygon(<ENTER VALUE HERE>) 'Amount of vertices in polygon - 1
For i = 0 To <ENTER VALUE HERE> 'Same value as above
Set polygon(i) = New CPoint
polygon(i).X = <ASSIGN X VALUE HERE> 'Source a list of values that can be looped through
polgyon(i).Y = <ASSIGN Y VALUE HERE> 'Source a list of values that can be looped through
Next i
InPolygon = isPointInPolygon(p, polygon)
MsgBox InPolygon
End Sub
入力
- bounding_points：ポリゴンを構成するノード。
-
bounding_box_positions：フィルタリングする候補ポイント。 （バウンディングボックスから作成された実装では。
（入力は次の形式のタプルのリストです：
[(xcord, ycord), ...]）
返品
- ポリゴンの内側にあるすべてのポイント。
def polygon_ray_casting(self, bounding_points, bounding_box_positions):
# Arrays containing the x- and y-coordinates of the polygon's vertices.
vertx = [point[0] for point in bounding_points]
verty = [point[1] for point in bounding_points]
# Number of vertices in the polygon
nvert = len(bounding_points)
# Points that are inside
points_inside = []
# For every candidate position within the bounding box
for idx, pos in enumerate(bounding_box_positions):
testx, testy = (pos[0], pos[1])
c = 0
for i in range(0, nvert):
j = i - 1 if i != 0 else nvert - 1
if( ((verty[i] > testy ) != (verty[j] > testy)) and
(testx < (vertx[j] - vertx[i]) * (testy - verty[i]) / (verty[j] - verty[i]) + vertx[i]) ):
c += 1
# If odd, that means that we are inside the polygon
if c % 2 == 1:
points_inside.append(pos)
return points_inside
繰り返しますが、アイデアはこちら
これは、レイキャスティングを使用しないCのポリゴンテストのポイントです。また、重複する領域（自己交差）でも機能します。use_holes引数を参照してください。
/* math lib (defined below) */
static float dot_v2v2(const float a[2], const float b[2]);
static float angle_signed_v2v2(const float v1[2], const float v2[2]);
static void copy_v2_v2(float r[2], const float a[2]);
/* intersection function */
bool isect_point_poly_v2(const float pt[2], const float verts[][2], const unsigned int nr,
const bool use_holes)
{
/* we do the angle rule, define that all added angles should be about zero or (2 * PI) */
float angletot = 0.0;
float fp1[2], fp2[2];
unsigned int i;
const float *p1, *p2;
p1 = verts[nr - 1];
/* first vector */
fp1[0] = p1[0] - pt[0];
fp1[1] = p1[1] - pt[1];
for (i = 0; i < nr; i++) {
p2 = verts[i];
/* second vector */
fp2[0] = p2[0] - pt[0];
fp2[1] = p2[1] - pt[1];
/* dot and angle and cross */
angletot += angle_signed_v2v2(fp1, fp2);
/* circulate */
copy_v2_v2(fp1, fp2);
p1 = p2;
}
angletot = fabsf(angletot);
if (use_holes) {
const float nested = floorf((angletot / (float)(M_PI * 2.0)) + 0.00001f);
angletot -= nested * (float)(M_PI * 2.0);
return (angletot > 4.0f) != ((int)nested % 2);
}
else {
return (angletot > 4.0f);
}
}
/* math lib */
static float dot_v2v2(const float a[2], const float b[2])
{
return a[0] * b[0] + a[1] * b[1];
}
static float angle_signed_v2v2(const float v1[2], const float v2[2])
{
const float perp_dot = (v1[1] * v2[0]) - (v1[0] * v2[1]);
return atan2f(perp_dot, dot_v2v2(v1, v2));
}
static void copy_v2_v2(float r[2], const float a[2])
{
r[0] = a[0];
r[1] = a[1];
}
注：atan2fへの呼び出しが多く含まれているため、これは最適ではないメソッドの1つですが、このスレッドを読んでいる開発者には興味深いかもしれません（私のテストでは、ライン交差メソッドを使用した場合〜23倍遅くなります） ）。
ポリゴンのヒットを検出するには、2つのことをテストする必要があります：
- ポイントがポリゴン領域内にある場合。 （レイキャスティングアルゴリズムで実現可能）
- ポイントがポリゴンの境界線上にある場合（ポリライン（ライン）上のポイント検出に使用されるのと同じアルゴリズムで実現できます）。
の次の特殊なケースに対処するにはレイキャスティングアルゴリズム：
- レイはポリゴンの側面の1つと重なります。
- ポイントはポリゴンの内側にあり、光線はポリゴンの頂点を通過します。
- ポイントはポリゴンの外側にあり、光線はポリゴンの角度の1つに接触します。
ポイントが複雑なポリゴン内にあるかどうかを確認します。この記事では、これらの問題を解決する簡単な方法を提供するため、上記の場合に特別な処理は必要ありません。
これを行うには、目的のポイントをポリゴンの頂点に接続して形成された領域が、ポリゴン自体の領域と一致するかどうかを確認します。
または、あなたのポイントからチェックポイントまでの2つの連続したポリゴン頂点の各ペアまでの内角の合計が360になるかどうかを確認することもできますが、最初のオプションは含まれていないため、最初のオプションが速いと感じています除算も三角関数の逆関数の計算も。
ポリゴンに内部に穴がある場合、どうなるかわかりませんが、主な考えはこの状況に適応できるようです
数学コミュニティで質問を投稿することもできます。 100万通りの方法があると思う
java-scriptライブラリを探している場合、Polygonクラス用のJavaScript googleマップv3拡張があり、その中にポイントが存在するかどうかを検出します。
var polygon = new google.maps.Polygon([], "#000000", 1, 1, "#336699", 0.3);
var isWithinPolygon = polygon.containsLatLng(40, -90);
qt を使用する場合（Qt 4.3以降） 、QPolygonの関数 containsPoint
を使用できます。答えは、単純なポリゴンか複雑なポリゴンかによって異なります。単純なポリゴンには、ラインセグメントの交差があってはなりません。そのため、穴を開けることはできますが、線は互いに交差できません。複雑な領域には、線の交点を含めることができます。そのため、重なり合う領域、または1点だけで互いに接触する領域を持つことができます。
単純なポリゴンの場合、最適なアルゴリズムはレイキャスティング（交差数）アルゴリズムです。複雑なポリゴンの場合、このアルゴリズムは重複領域内にあるポイントを検出しません。したがって、複雑なポリゴンの場合は、巻数アルゴリズムを使用する必要があります。
これは、両方のアルゴリズムのC実装に関する優れた記事です。私はそれらを試してみましたが、うまくいきます。
nirgによるソリューションのスカラバージョン（境界矩形の事前チェックが個別に行われると想定）：
def inside(p: Point, polygon: Array[Point], bounds: Bounds): Boolean = {
val length = polygon.length
@tailrec
def oddIntersections(i: Int, j: Int, tracker: Boolean): Boolean = {
if (i == length)
tracker
else {
val intersects = (polygon(i).y > p.y) != (polygon(j).y > p.y) && p.x < (polygon(j).x - polygon(i).x) * (p.y - polygon(i).y) / (polygon(j).y - polygon(i).y) + polygon(i).x
oddIntersections(i + 1, i, if (intersects) !tracker else tracker)
}
}
oddIntersections(0, length - 1, tracker = false)
}
@nirg answerのgolangバージョンです（@@ m-katzのC＃コードに触発された）
func isPointInPolygon(polygon []point, testp point) bool {
minX := polygon[0].X
maxX := polygon[0].X
minY := polygon[0].Y
maxY := polygon[0].Y
for _, p := range polygon {
minX = min(p.X, minX)
maxX = max(p.X, maxX)
minY = min(p.Y, minY)
maxY = max(p.Y, maxY)
}
if testp.X < minX || testp.X > maxX || testp.Y < minY || testp.Y > maxY {
return false
}
inside := false
j := len(polygon) - 1
for i := 0; i < len(polygon); i++ {
if (polygon[i].Y > testp.Y) != (polygon[j].Y > testp.Y) && testp.X < (polygon[j].X-polygon[i].X)*(testp.Y-polygon[i].Y)/(polygon[j].Y-polygon[i].Y)+polygon[i].X {
inside = !inside
}
j = i
}
return inside
}
これを早めに提起した人はいませんでしたが、データベースを必要とするプラグマティスト向け：MongoDBは、これを含むGeoクエリに対して優れたサポートを提供しています。
探しているのは：
db.neighborhoods.findOne（{geometry：{$ geoIntersects：{$ geometry：{ タイプ：<！> quot;ポイント<！> quot ;、座標：[<！> quot; longitude <！> quot;、<！> quot; latitude <！> quot; ]}}} }）
Neighborhoodsは、1つ以上のポリゴンを標準のGeoJson形式で保存するコレクションです。クエリがnullを返した場合は交差しません。そうでない場合は交差します。
非常によくここに文書化されています： https://docs.mongodb.com/manual/tutorial/geospatial-tutorial/
330の不規則なポリゴングリッドに分類された6,000ポイントを超えるパフォーマンスは、最適化がまったく行われず、ドキュメントをそれぞれのポリゴンで更新する時間も含めて1分未満でした。
