我试图创建一个快速的2D点内多边形算法,用于命中测试(例如多边形.contains(p:点))。对有效技术的建议将不胜感激。
当前回答
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;
}
其他回答
简单的解决方案是将多边形划分为三角形,并按这里解释的那样对三角形进行测试
如果你的多边形是凸多边形,可能有更好的方法。把这个多边形看作是无限条线的集合。每一行将空间一分为二。对于每一个点,很容易判断它是在直线的一边还是另一边。如果一个点在所有直线的同一侧,那么它在多边形内。
当我还是Michael Stonebraker手下的一名研究员时,我做了一些关于这方面的工作——你知道,就是那位提出了Ingres、PostgreSQL等的教授。
我们意识到最快的方法是首先做一个边界框,因为它非常快。如果它在边界框之外,它就在外面。否则,你就得做更辛苦的工作……
如果你想要一个伟大的算法,看看开源项目PostgreSQL的源代码的地理工作…
我想指出的是,我们从来没有深入了解过左撇子和右撇子(也可以表达为“内”和“外”的问题……
更新
BKB's link provided a good number of reasonable algorithms. I was working on Earth Science problems and therefore needed a solution that works in latitude/longitude, and it has the peculiar problem of handedness - is the area inside the smaller area or the bigger area? The answer is that the "direction" of the verticies matters - it's either left-handed or right handed and in this way you can indicate either area as "inside" any given polygon. As such, my work used solution three enumerated on that page.
此外,我的工作使用单独的函数进行“在线”测试。
...因为有人问:我们发现当垂直的数量超过某个数字时,边界盒测试是最好的——如果有必要,在做更长的测试之前做一个非常快速的测试……边界框是通过简单地将最大的x,最小的x,最大的y和最小的y放在一起,组成一个框的四个点来创建的……
另一个提示是:我们在网格空间中进行了所有更复杂的“调光”计算,都是在平面上的正点上进行的,然后重新投影到“真实”的经度/纬度上,从而避免了在经度180线交叉时和处理极地时可能出现的环绕错误。工作好了!
下面是Rust版本的@nirg答案(Philipp Lenssen javascript版本) 我给出这个答案是因为我从这个网站得到了很多帮助,我翻译javascript版本rust作为一个练习,希望可以帮助一些人,最后一个原因是,在我的工作中,我会把这段代码翻译成一个wasm,以提高我的画布的性能,这是一个开始。我的英语很差……,请原谅我 `
pub struct Point {
x: f32,
y: f32,
}
pub fn point_is_in_poly(pt: Point, polygon: &Vec<Point>) -> bool {
let mut is_inside = false;
let max_x = polygon.iter().map(|pt| pt.x).reduce(f32::max).unwrap();
let min_x = polygon.iter().map(|pt| pt.x).reduce(f32::min).unwrap();
let max_y = polygon.iter().map(|pt| pt.y).reduce(f32::max).unwrap();
let min_y = polygon.iter().map(|pt| pt.y).reduce(f32::min).unwrap();
if pt.x < min_x || pt.x > max_x || pt.y < min_y || pt.y > max_y {
return is_inside;
}
let len = polygon.len();
let mut j = len - 1;
for i in 0..len {
let y_i_value = polygon[i].y > pt.y;
let y_j_value = polygon[j].y > pt.y;
let last_check = (polygon[j].x - polygon[i].x) * (pt.y - polygon[i].y)
/ (polygon[j].y - polygon[i].y)
+ polygon[i].x;
if y_i_value != y_j_value && pt.x < last_check {
is_inside = !is_inside;
}
j = i;
}
is_inside
}
let pt = Point {
x: 1266.753,
y: 97.655,
};
let polygon = vec![
Point {
x: 725.278,
y: 203.586,
},
Point {
x: 486.831,
y: 441.931,
},
Point {
x: 905.77,
y: 445.241,
},
Point {
x: 1026.649,
y: 201.931,
},
];
let pt1 = Point {
x: 725.278,
y: 203.586,
};
let pt2 = Point {
x: 872.652,
y: 321.103,
};
println!("{}", point_is_in_poly(pt, &polygon));// false
println!("{}", point_is_in_poly(pt1, &polygon)); // true
println!("{}", point_is_in_poly(pt2, &polygon));// true
`
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;
}
}
没有什么比归纳定义问题更美好的了。为了完整起见,你在序言中有一个版本,它可能也澄清了光线投射背后的思想:
基于仿真的简化算法在http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
一些helper谓词:
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的直线(直线(a,B))方程为:
(YB-YA)
Y - YA = ------- * (X - XA)
(XB-YB)
It is important that the direction of rotation for the line is setted to clock-wise for boundaries and anti-clock-wise for holes. We are going to check whether the point (X,Y), i.e the tested point is at the left half-plane of our line (it is a matter of taste, it could also be the right side, but also the direction of boundaries lines has to be changed in that case), this is to project the ray from the point to the right (or left) and acknowledge the intersection with the line. We have chosen to project the ray in the horizontal direction (again it is a matter of taste, it could also be done in vertical with similar restrictions), so we have:
(XB-XA)
X < ------- * (Y - YA) + XA
(YB-YA)
Now we need to know if the point is at the left (or right) side of the line segment only, not the entire plane, so we need to restrict the search only to this segment, but this is easy since to be inside the segment only one point in the line can be higher than Y in the vertical axis. As this is a stronger restriction it needs to be the first to check, so we take first only those lines meeting this requirement and then check its possition. By the Jordan Curve theorem any ray projected to a polygon must intersect at an even number of lines. So we are done, we will throw the ray to the right and then everytime it intersects a line, toggle its state. However in our implementation we are goint to check the lenght of the bag of solutions meeting the given restrictions and decide the innership upon it. for each line in the polygon this have to be done.
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).
推荐文章
- 如何加速gwt编译器?
- MySQL OR与IN性能
- 应该……接住环内还是环外?
- 哪个更快/最好?SELECT *或SELECT columnn1, colum2, column3等
- 加快R中的循环操作
- INT和VARCHAR主键之间有真正的性能差异吗?
- c++标准是否要求iostreams的性能很差,或者我只是在处理一个糟糕的实现?
- 大概的成本访问各种缓存和主存储器?
- 模拟慢速互联网连接
- 如何检查表上持有哪些锁
- 检查字符串是否包含字符串列表中的元素
- Scala vs Python的Spark性能
- 现代c++能让你免费获得性能吗?
- 对于PostgreSQL表来说,多大才算太大?
- 即使从未抛出异常,使用try-catch块的代价是否昂贵?