ヒットテスト (例) で使用するために、ポリゴン アルゴリズム内に高速な 2D ポイントを作成しようとしていますPolygon.contains(p:Point)
。効果的なテクニックの提案をいただければ幸いです。
38 に答える
グラフィックの場合、私はむしろ整数を好みません。多くのシステムは UI の描画に整数を使用しますが (結局ピクセルは int です)、たとえば macOS はすべてに float を使用します。macOS はポイントのみを認識し、ポイントは 1 ピクセルに変換できますが、モニターの解像度によっては別のピクセルに変換される場合があります。Retina スクリーンでは、ポイントの半分 (0.5/0.5) がピクセルです。それでも、macOS の UI が他の UI よりも大幅に遅いことに気づいたことはありません。結局、3D API (OpenGL または Direct3D) は浮動小数点数でも動作し、最新のグラフィックス ライブラリは GPU アクセラレーションを頻繁に利用します。
速度が主な関心事であると言いましたが、速度を上げましょう。高度なアルゴリズムを実行する前に、まず簡単なテストを行います。ポリゴンの周りに軸に沿ったバウンディング ボックスを作成します。これは非常に簡単で高速であり、すでに多くの計算を節約できます。それはどのように機能しますか?多角形のすべてのポイントを反復処理し、X と Y の最小値/最大値を見つけます。
たとえば、あなたはポイントを持っています(9/1), (4/3), (2/7), (8/2), (3/6)
。これは Xmin が 2、Xmax が 9、Ymin が 1、Ymax が 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 回だけ当たるか、まったく当たらない必要があります。同じ面に 2 回ヒットすることはありません。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;
}
入力値は、ベクトル 1 (および) とベクトル 2 (および) の2 つのエンドポイントです。つまり、2 つのベクトル、4 つの点、8 つの座標があります。そして明確です。交差を増やしますが、何もしません。v1x1/v1y1
v1x2/v1y2
v2x1/v2y1
v2x2/v2y2
YES
NO
YES
NO
COLLINEARはどうですか?これは、位置と長さに応じて、両方のベクトルが同じ無限の線上にあることを意味し、まったく交差しないか、無限の数の点で交差します。このケースをどのように処理するかは完全にはわかりません。どちらの方法でも交差点としてカウントしません。まあ、浮動小数点の丸め誤差のため、このケースは実際にはかなりまれです。より良いコードはおそらくテストでは== 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座標。
短くて効率的で、凸多角形と凹多角形の両方で機能します。前に提案したように、最初に外接する長方形を確認し、ポリゴンの穴を個別に処理する必要があります。
The idea behind this is pretty simple. The author describes it as follows:
I run a semi-infinite ray horizontally (increasing x, fixed y) out from the test point, and count how many edges it crosses. At each crossing, the ray switches between inside and outside. This is called the Jordan curve theorem.
The variable c is switching from 0 to 1 and 1 to 0 each time the horizontal ray crosses any edge. So basically it's keeping track of whether the number of edges crossed are even or odd. 0 means even and 1 means odd.
これは、この RPI 教授からのnirg によって与えられた回答の C# バージョンです。その 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;
}
// https://wrf.ecse.rpi.edu/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;
}
点 p と各多角形の頂点の間の角度の合計を計算します。方向角の合計が 360 度の場合、ポイントは内側にあります。合計が 0 の場合、ポイントは外側です。
この方法の方が堅牢で、数値精度への依存度が低いため、私はこの方法の方が気に入っています。
交差数の計算中に頂点に「ヒット」する可能性があるため、交差数の均一性を計算する方法は限られています。
EDIT:ちなみに、この方法は凹面ポリゴンと凸面ポリゴンで機能します。
編集: 最近、このトピックに関するウィキペディアの記事全体を見つけました。
ボボボボが引用したエリック・ヘインズの記事は本当に素晴らしいです。特に興味深いのは、アルゴリズムのパフォーマンスを比較した表です。角度の合計方法は、他の方法に比べて本当に悪いです。また興味深いのは、ルックアップ グリッドを使用してポリゴンをさらに "in" セクターと "out" セクターに分割するなどの最適化により、1000 面を超えるポリゴンでもテストが非常に高速になることです。
とにかく、まだ始まったばかりですが、私の投票は「クロッシング」方法に行きます。しかし、David Bourke によって最も簡潔に記述され体系化されていることがわかりました。実際の三角法は必要なく、凸面でも凹面でも機能し、辺の数が増えるにつれて適度に機能することが気に入っています。
ちなみに、Eric Haines の記事のパフォーマンス テーブルの 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
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の下で研究者だったときに、これについていくつかの作業を行いました。
超高速なので、最初に境界ボックスを作成するのが最速の方法であることに気付きました。境界ボックスの外側にある場合は、外側です。そうしないと、あなたはもっと難しい仕事をします...
優れたアルゴリズムが必要な場合は、オープン ソース プロジェクトの PostgreSQL ソース コードを参照して、地理的な作業を行ってください...
指摘したいのは、右利きと左利きについての洞察はまったく得られていないことです(「内側」と「外側」の問題としても表現できます...
アップデート
BKB のリンクは、妥当なアルゴリズムを数多く提供しています。私は地球科学の問題に取り組んでいたため、緯度/経度で機能するソリューションが必要でした。それには利き手という特有の問題があります-小さいエリアまたは大きいエリアの内側のエリアですか? 答えは、頂点の「方向」が重要であるということです。左手または右手のいずれかであり、このようにして、特定のポリゴンの「内側」としていずれかの領域を示すことができます。そのため、私の作業では、そのページに列挙されているソリューション 3 を使用しました。
さらに、私の作品では、「オンライン」テスト用に別の関数を使用していました。
...誰かが尋ねたので: 頂点の数がある数を超えたときにバウンディング ボックス テストが最適であることがわかりました.必要に応じて長いテストを行う前に非常に簡単なテストを行ってください.最大の x、最小の x、最大の y、最小の y を組み合わせて、ボックスの 4 つのポイントを作成します...
続く人のためのもう 1 つのヒント: グリッド空間ですべてのより洗練された「ライトディミング」コンピューティングをすべて平面上の正の点で実行し、「実際の」経度/緯度に再投影して、エラーの可能性を回避しました。経度 180 度の線を横切ったとき、および極地域を処理するときにラップアラウンドします。うまくいきました!
David Segond の回答は標準的な一般的な回答であり、Richard T の回答は最も一般的な最適化ですが、他にもいくつかあります。他の強力な最適化は、一般的ではないソリューションに基づいています。たとえば、多数のポイントを持つ同じポリゴンをチェックする場合、非常に高速な TIN 検索アルゴリズムが多数あるため、ポリゴンを三角測量すると処理速度が大幅に向上します。もう1つは、ポリゴンとポイントが低解像度の限られた平面上にある場合、たとえば画面表示の場合、ポリゴンを特定の色でメモリマップされたディスプレイバッファーにペイントし、特定のピクセルの色をチェックして、それが存在するかどうかを確認できますポリゴンで。
多くの最適化と同様に、これらは一般的なケースではなく特定のケースに基づいており、1 回の使用ではなく償却された時間に基づいてメリットが得られます。
この分野で働いていると、Jooseph O'Rourkes の「Computation Geometry in C」ISBN 0-521-44034-3 が大きな助けになることがわかりました。
簡単な解決策は、ポリゴンを三角形に分割し、ここで説明されているように三角形のヒット テストを行うことです。
ポリゴンが凸面の場合は、より良いアプローチがあるかもしれません。多角形を無限の線の集まりとして見てください。スペースを 2 つに分割する各線。すべての点について、それが線の片側にあるか反対側にあるかを簡単に言うことができます。ポイントがすべてのラインの同じ側にある場合、そのポイントはポリゴンの内側にあります。
これは古いことだと思いますが、誰かが興味を持っている場合に備えて、Cocoaに実装されているレイキャスティングアルゴリズムを次に示します。それが物事を行うための最も効率的な方法であるかどうかはわかりませんが、誰かを助けるかもしれません。
- (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;
}
問題の帰納的な定義ほど美しいものはありません。ここで完全を期すために、プロローグにバージョンがあり、レイキャスティングの背後にある考えも明確にする可能性があります。
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]).
2 つの点 A と B が与えられた直線の方程式 (Line(A,B)) は次のとおりです。
(YB-YA)
Y - YA = ------- * (X - XA)
(XB-YB)
ラインの回転方向は、境界では時計回り、穴では反時計回りに設定することが重要です。点 (X,Y)、つまりテストされた点が線の左半平面にあるかどうかを確認します (好みの問題です。右側である可能性もありますが、境界の方向でもある可能性があります)。その場合は線を変更する必要があります)、これは点から右 (または左) に光線を投影し、線との交点を認識するためです。光線を水平方向に投影することを選択したので (これも好みの問題ですが、同様の制限で垂直方向にも行うことができます)、次のようになります。
(XB-XA)
X < ------- * (Y - YA) + XA
(YB-YA)
ここで、点が平面全体ではなく、線分の左側 (または右側) にあるかどうかを知る必要があるため、検索をこの線分のみに制限する必要がありますが、線分内にあるため、これは簡単です。縦軸の Y よりも高い位置にできるのは、ライン内の 1 つのポイントだけです。これはより強力な制限であるため、最初にチェックする必要があるため、最初にこの要件を満たす行のみを取得し、その位置を確認します。ジョーダン曲線の定理により、多角形に射影された光線は、偶数の線で交差する必要があります。これで完了です。光線を右に投げ、線と交差するたびにその状態を切り替えます。ただし、実装では、指定された制限を満たすソリューションのバッグの長さをチェックし、それに対する内部関係を決定します。ポリゴンの各ラインに対して、これを行う必要があります。
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).
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;
}
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
.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;
}
}
レイキャスティングを使用していない 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];
}
注: これは への呼び出しが多く含まれているため、あまり最適ではない方法の 1 つですが、atan2f
このスレッドを読んでいる開発者にとっては興味深いかもしれません (私のテストでは、線交差法を使用する場合よりも 23 倍遅くなります)。
これはおそらく、このページから入手した hereの C コードのわずかに最適化されていないバージョンです。
私の C++ バージョンではstd::vector<std::pair<double, double>>
、x と y として a と 2 つの double を使用します。ロジックは元の C コードとまったく同じはずですが、私のほうが読みやすいと思います。パフォーマンスについて話すことはできません。
bool point_in_poly(std::vector<std::pair<double, double>>& verts, double point_x, double point_y)
{
bool in_poly = false;
auto num_verts = verts.size();
for (int i = 0, j = num_verts - 1; i < num_verts; j = i++) {
double x1 = verts[i].first;
double y1 = verts[i].second;
double x2 = verts[j].first;
double y2 = verts[j].second;
if (((y1 > point_y) != (y2 > point_y)) &&
(point_x < (x2 - x1) * (point_y - y1) / (y2 - y1) + x1))
in_poly = !in_poly;
}
return in_poly;
}
元のCコードは
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;
}
レイ キャスティング アルゴリズムで次の特殊なケースに対処するには:
- レイは、ポリゴンの側面の 1 つに重なっています。
- 点は多角形の内側にあり、光線は多角形の頂点を通過します。
- ポイントはポリゴンの外側にあり、レイはポリゴンの角度の 1 つにちょうど接触します。
点が複雑な多角形の内側にあるかどうかの判断を確認します。この記事では、それらを解決する簡単な方法を提供しているため、上記のケースに特別な処置は必要ありません。
java-script ライブラリを探している場合は、Polygon クラスの javascript google maps v3 拡張機能があり、ポイントがその中に存在するかどうかを検出します。
var polygon = new google.maps.Polygon([], "#000000", 1, 1, "#336699", 0.3);
var isWithinPolygon = polygon.containsLatLng(40, -90);
これは凸形状に対してのみ機能しますが、ミンコフスキー ポータル リファインメントと GJK は、ポイントがポリゴン内にあるかどうかをテストするための優れたオプションでもあります。ミンコフスキー減算を使用してポリゴンからポイントを減算し、それらのアルゴリズムを実行して、ポリゴンに原点が含まれているかどうかを確認します。
また、興味深いことに、方向ベクトルを入力として取り、そのベクトルに沿って最も遠い点を吐き出すサポート関数を使用して、形状をもう少し暗黙的に記述することができます。これにより、曲線、多角形、混合など、あらゆる凸形状を表現できます。単純なサポート関数の結果を組み合わせて、より複雑な形状を作成する操作を実行することもできます。
詳細: http://xenocollide.snethen.com/mpr2d.html
また、Game Programming gems 7 では、これを 3D で行う方法について説明しています (: