如何确定2D点是否在多边形内？
问题描述 投票：458回答：32
我正在尝试在多边形算法中创建一个快速的2D点，用于命中测试（例如Polygon.contains(p:Point)）。对于有效技术的建议将不胜感激。
32个回答
投票
对于图形，我宁愿不喜欢整数。许多系统使用整数进行UI绘制（毕竟像素是整数），但macOS例如使用float来表示一切。 macOS只知道点，一个点可以转换为一个像素，但根据显示器的分辨率，它可能会转换为其他像素。在视网膜屏幕上，半个点（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.具有两个边（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!
}
这是第一个针对任何一点运行的测试。正如您所看到的，此测试速度非常快，但也非常粗糙。要处理边界矩形内的点，我们需要更复杂的算法。有几种方法可以计算出来。如果多边形可以有孔或者总是是实心的，那么哪种方法也起作用。以下是实体（一个凸面，一个凹面）的示例：
而这里有一个洞：
绿色的中间有一个洞！
最简单的算法，可以处理上述所有三种情况，并且仍然非常快，称为光线投射。该算法的想法非常简单：从多边形外部的任何位置绘制虚拟光线到您的点，并计算它撞击多边形一侧的频率。如果命中数是偶数，则它在多边形之外，如果它是奇数，则它在内部。
绕组数算法是一种替代方案，它对于非常接近多边形线的点更准确，但它也慢得多。由于有限的浮点精度和圆角问题，射线投射可能因太靠近多边形的点而失败，但实际上这几乎不是问题，就好像一个点靠近一侧，它通常在视觉上甚至不可能观众识别它是否已经在里面或者仍然在外面。
你还有上面的边框，还记得吗？只需在边界框外选取一个点，并将其用作光线的起点。例如。点(Xmin - e/p.y)肯定在多边形之外。
但e是什么？好吧，e（实际上是epsilon）为边界框提供了一些填充。正如我所说，如果我们太靠近多边形线，光线跟踪就会失败。由于边界框可能等于多边形（如果多边形是一个轴对齐的矩形，边界框等于多边形本身！），我们需要一些填充来使这个安全，就是这样。你应该选择e多大？不太大。它取决于您用于绘图的坐标系比例。如果您的像素步长为1.0，那么只需选择1.0（但是0.1也会有效）
现在我们有了射线的起点和终点坐标，问题从“多边形内的点”变为“射线与多边形面相交的频率”。因此，我们不能像以前一样使用多边形点，现在我们需要实际的边。一方总是由两点来定义。
side 1: (X1/Y1)-(X2/Y2)
side 2: (X2/Y2)-(X3/Y3)
side 3: (X3/Y3)-(X4/Y4)
:
你需要测试所有侧面的光线。将光线视为矢量，将每一边视为矢量。射线必须完全击中每一侧或从不击中每一侧。它不能两次击中同一侧。 2D空间中的两条线总是恰好相交一次，除非它们是平行的，在这种情况下它们永远不会相交。然而，由于矢量具有有限的长度，因此两个矢量可能不是平行的并且仍然不会相交，因为它们太短而不能彼此相遇。
// 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
}
到目前为止一切顺利，但是你如何测试两个向量是否相交？这里有一些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（v1x1/v1y1和v1x2/v1y2）和向量2（v2x1/v2y1和v2x2/v2y2）的两个端点。所以你有2个向量，4个点，8个坐标。 YES和NO很清楚。 YES增加交叉点，NO什么都不做。
COLLINEAR怎么样？这意味着两个矢量位于相同的无限线上，取决于位置和长度，它们根本不相交或者它们在无数个点中相交。我不完全确定如何处理这种情况，我不认为它是交叉点。嗯，这种情况在实践中相当罕见，因为浮点舍入错误;更好的代码可能不会测试== 0.0f，而是像< epsilon，epsilon是一个相当小的数字。
如果你需要测试更多的点，你可以通过将多边形边的线性方程式标准形式保存在内存中来加快整个过程，所以你不必每次都重新计算这些。这将在每次测试中为您节省两次浮点乘法和三次浮点减法，以换取在内存中存储每个多边形侧的三个浮点值。这是典型的内存与计算时间的权衡。
最后但并非最不重要：如果您可以使用3D硬件来解决问题，那么有一个有趣的选择。让GPU为您完成所有工作。创建一个屏幕外的绘画表面。用黑色完全填充。现在让OpenGL或Direct3D绘制多边形（或者甚至是所有多边形，如果你只想测试点是否在任何一个内，但你不关心哪一个）并用不同的方法填充多边形颜色，例如白色。要检查点是否在多边形内，请从绘图表面获取此点的颜色。这只是一个O（1）内存提取。
当然，只有在绘图表面不必很大的情况下，此方法才可用。如果它不能适合GPU内存，这种方法比在CPU上运行要慢。如果它必须是巨大的并且您的GPU支持现代着色器，您仍然可以通过将上面显示的光线投射实现为GPU着色器来使用GPU，这绝对是可能的。对于更大数量的多边形或大量的测试点，这将得到回报，考虑一些GPU将能够并行测试64到256个点。但请注意，将数据从CPU传输到GPU并返回总是很昂贵，因此只需针对几个简单的多边形测试几个点，其中点或多边形是动态的并且会经常变化，GPU方法很少会付出代价关闭。
投票
真的很喜欢由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搜索算法。另一个是如果多边形和点在低分辨率的有限平面上，比如屏幕显示，您可以将多边形绘制到给定颜色的内存映射显示缓冲区中，并检查给定像素的颜色以查看它是否位于在多边形中。
与许多优化一样，这些优化基于特定情况而非一般情况，并且基于摊销时间而非单次使用产生收益。
在这个领域工作，我发现Joseph O'Rourke的'计算几何在C'国际标准书号0-521-44034-3是一个很好的帮助。
投票
平凡的解决方案是将多边形划分为三角形并按照here的说明对三角形进行测试
如果您的多边形是CONVEX，那么可能有更好的方法。将多边形视为无限线的集合。每条线将空间分成两部分。对于每一点，很容易说它是在线的一侧还是另一侧。如果一个点位于所有线的同一侧，则它位于多边形内。
投票
我意识到这是旧的，但这是一个在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;
}
投票
Obir-C版本的nirg的答案与测试点的样本方法。 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;
}
投票
没有什么比问题的归纳定义更美观了。为了完整起见，你在prolog中有一个版本，它也可能澄清了光线投射背后的想法：
基于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（线（A，B））的线的方程是：
(YB-YA)
Y - YA = ------- * (X - XA)
(XB-YB)
重要的是，线的旋转方向设置为时钟方向的边界和反时钟方式的孔。我们要检查点（X，Y），即测试点是否在我们线的左半平面（这是一个品味的问题，它也可能是右侧，但也是边界的方向在这种情况下必须改变线条），这是将光线从点投射到右边（或左边）并确认与线的交点。我们选择在水平方向上投射光线（这也是一个品味问题，它也可以在垂直方向上进行，具有类似的限制），因此我们有：
(XB-XA)
X < ------- * (Y - YA) + XA
(YB-YA)
现在我们需要知道点是仅位于线段的左侧（或右侧），而不是整个平面，因此我们需要将搜索仅限制到此段，但这很容易，因为它位于段内在垂直轴上，线中只有一个点可以高于Y.由于这是一个更强的限制，它需要首先检查，所以我们首先只采取符合此要求的那些线，然后检查它的位置。根据Jordan Curve定理，投射到多边形的任何光线必须在偶数行处相交。所以我们完成了，我们将光线投射到右边，然后每次它与一条线相交，切换它的状态。然而，在我们的实施中，我们要检查满足给定限制的解决方案包的长度，并决定其中的内容。对于多边形中的每一行，必须这样做。
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版本：
注意：请记住，如果您的多边形是地图中的一个区域，纬度/经度是Y / X值而不是X / Y（纬度= Y，经度= X），因为据我所知，这是历史的影响。经度不是衡量标准。
CLASS MODULE：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
投票
我认为以下代码是最好的解决方案（取自here）：
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;
}
Arguments
- convert：多边形中的顶点数。在上面提到的文章中已经讨论了是否在末尾重复第一个顶点。
- vertx，verty：包含多边形顶点的x坐标和y坐标的数组。
- testx，testy：测试点的X坐标和y坐标。
它既短又高效，适用于凸多边形和凹多边形。如前所述，您应首先检查边界矩形并分别处理多边形孔。
这背后的想法非常简单。作者将其描述如下：
我从测试点水平运行半无限光线（增加x，固定y），并计算它穿过的边数。在每个交叉点，射线在内部和外部之间切换。这被称为乔丹曲线定理。
每当水平光线穿过任何边缘时，变量c从0切换到1并且从1切换到0。所以基本上它是跟踪交叉的边数是偶数还是奇数。 0表示偶数，1表示奇数。
投票
我已经制作了nirg's c ++ code的Python实现：
输入
- 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
同样，这个想法取自here
投票
在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的调用，但是读取这个线程的开发人员可能会感兴趣（在我的测试中，它比使用行交集方法慢23倍）。
投票
对于检测多边形上的命中，我们需要测试两件事：
- 如果Point在多边形区域内。 （可以通过射线投射算法完成）
- 如果Point位于多边形边界上（可以通过用于折线（线）上的点检测的相同算法完成）。
投票
处理Ray casting algorithm中的以下特殊情况：
- 光线与多边形的一侧重叠。
- 该点位于多边形内部，并且光线穿过多边形的顶点。
- 该点位于多边形之外，并且光线刚好接触多边形的一个角度。
检查Determining Whether A Point Is Inside A Complex Polygon。该文章提供了一种解决问题的简便方法，因此上述案例不需要特殊处理。
投票
您可以通过检查通过将所需点连接到多边形的顶点而形成的区域是否与多边形本身的区域匹配来执行此操作。
或者您可以检查从您的点到每对两个连续多边形顶点到检查点的内角之和是否总和为360，但我感觉第一个选项更快，因为它不涉及分割或计算三角函数的逆。
我不知道如果你的多边形里面有一个洞会发生什么，但在我看来，主要的想法可以适应这种情况
您也可以在数学社区中发布问题。我打赌他们有一百万种方法可以做到这一点
投票
如果您正在寻找一个java脚本库，那么Polygon类的javascript google maps 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
投票
答案取决于你是否有简单或复杂的多边形。简单多边形不得有任何线段交叉点。所以他们可以有洞，但线不能相互交叉。复杂区域可以具有线交叉 - 因此它们可以具有重叠区域，或者仅通过单个点彼此接触的区域。
对于简单多边形，最好的算法是Ray cast（交叉数）算法。对于复杂多边形，此算法不会检测重叠区域内的点。因此，对于复杂多边形，您必须使用绕组数算法。
这是两篇算法的C实现的优秀文章。我尝试过它们并且效果很好。
投票
nirg解决方案的Scala版本（假定边界矩形预检是单独完成的）：
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)
}
投票
这是golang版的@nirg答案（灵感来自@@ 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
}
投票
这是answer given by nirg的C＃版本，来自this RPI professor。请注意，使用该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;
}
投票
令人惊讶的是没有人提前提出这个问题，但对于需要数据库的实用主义者来说：MongoDB对Geo查询提供了极好的支持，包括这个。
你在寻找的是：
db.neighborhoods.findOne（{geometry：{$ geoIntersects：{$ geometry：{type：“Point”，coordinates：[“longitude”，“latitude”]}}}}}）
Neighborhoods是以标准GeoJson格式存储一个或多个多边形的集合。如果查询返回null，则它不会相交，否则就是。
这里有很好的记录：https://docs.mongodb.com/manual/tutorial/geospatial-tutorial/
分类在330不规则多边形网格中的超过6,000个点的性能不到一分钟，根本没有优化，包括用各自的多边形更新文档的时间。
投票
以下是M. Katz基于Nirg方法的答案的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，则该点在外面。
我更喜欢这种方法，因为它更稳健，更少依赖于数值精度。
计算交叉点数均匀度的方法是有限的，因为您可以在计算交叉点数量时“击中”顶点。
编辑：通过The Way，这种方法适用于凹凸多边形。
编辑：我最近在这个主题上发现了一个完整的Wikipedia article。
投票
bobobobo引用的Eric Haines article非常出色。特别有趣的是比较算法性能的表格;角度求和方法与其他方法相比非常糟糕。同样有趣的是，使用查找网格进一步将多边形细分为“in”和“out”扇区的优化可以使测试速度极快，即使在具有> 1000边的多边形上也是如此。
无论如何，这是早期的，但我的投票是针对“交叉”方法，这正是Mecki所描述的。然而，我发现它最无助的described and codified by David Bourke。我喜欢不需要真正的三角法，它适用于凸面和凹面，并且随着边数的增加，它表现得相当好。
顺便说一句，这是Eric Haines的文章中的一个性能表，用于测试随机多边形。
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
投票
这个问题非常有趣。我有另一个可行的想法，不同于这篇文章的其他答案。
我们的想法是使用角度之和来确定目标是在内部还是外部。如果目标位于区域内，则目标和每两个边界点形成的角度之和将为360.如果目标位于外部，则总和将不是360.角度具有方向。如果角度向后，角度是负角度。这就像计算winding number一样。
请参阅此图片以基本了解该想法：
我的算法假设顺时针是正方向。这是一个潜在的输入：
[[-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的地理工作源代码......
我想指出，我们从来没有对左右手的任何见解（也可以表达为“内部”与“外部”问题......
UPDATE
BKB的链接提供了大量合理的算法。我正在研究地球科学问题，因此需要一个在纬度/经度上工作的解决方案，它有一个特殊的手性问题 - 是较小区域内还是较大区域内的区域？答案是顶点的“方向”很重要 - 它既可以是左手也可以是右手，这样你就可以将任一区域指定为任何给定多边形的“内部”。因此，我的工作使用了该页面上列举的解决方案三。
另外，我的工作使用了“在线”测试的单独功能。
...自从有人问：我们发现，当顶点数量超过一定数量时，边界框测试最好 - 如果需要，在进行更长时间的测试之前做一个非常快速的测试...通过简单地获取边界框来创建边界框最大x，最小x，最大y和最小y并将它们组合在一起制作一个盒子的四个点......
对于接下来的人的另一个提示：我们在网格空间中进行了所有更复杂和“光线调暗”的计算，所有这些都在平面上的正点上，然后重新投射回“真实的”经度/纬度，从而避免了可能的错误。当经过一条经线180和处理极地区域时，环绕着。工作得很好！
最新问题
- PHP 中来自 CSV / TXT 文件的关联数组
- Pandas 数据框 select_dtypes 未选择预期的数据类型
- clang AST 用于单行多变量声明
- SOME/IP Service Discovery中的Eventgroup ID从何而来
- COUNTA 的替代品
- Laravel artisan 命令 php artisan make:导入未定义
- 如何使用JNLP将EKS连接到ALB后面的jenkins？
- 如何旋转画布上的特定图像层？
- B2C 自定义策略上传时验证错误
- Box 实例仅复制到最后一行
- 内容处置内联文件名不起作用
- Flutter flutter_local_notifications 包 - 可以检查在线状态吗？
- Microsoft Azure：条件访问：如何强制执行 intune MDM 但允许来自个人设备的 AVD RDP？
- 更改边框宽度时如何防止移位
- 当前类的API文档预览
- 快速 UI 中自定义形状的条纹图案
- 学习新工具/库的正确程序是什么？ [已关闭]
- 如何在linux中限制c代码的堆大小
- 如何在 Swift 中为 UIView 子类编写自定义 init？
- 我应该如何为我的 AWS 无服务器应用程序迁移（更新数据库架构）我的数据库