我试图创建一个快速的2D点内多边形算法,用于命中测试(例如多边形.contains(p:点))。对有效技术的建议将不胜感激。
当前回答
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;
}
}
其他回答
在大多数情况下,这是一个比其他算法都快的算法。
它又新又雅致。我们花费O(n * log(n))时间构建一个表,允许我们在O(log(n) + k)时间内测试多边形中的点。
与光线跟踪或角度不同,使用扫描光束表可以更快地对同一多边形进行多次检查。我们必须预先构建一个扫描束活动边表,这是大多数代码正在做的事情。
We calculate the scanbeam and the active edges for that position in the y-direction. We make a list of points sorted by their y-component and we iterate through this list, for two events. Start-Y and End-Y, we track the active edges as we process the list. We process the events in order and for each scanbeam we record the y-value of the event and the active edges at each event (events being start-y and end-y) but we only record these when our event-y is different than last time (so everything at the event point is processed before we mark it in our table).
我们得到我们的表格:
[] p6p5、p6p7 p6p5, p6p7, p2p3, p2p1 p6p7, p5p4, p2p3, p3p1 p7p8, p5p4, p2p3, p2p1 p7p8, p5p4, p3p4, p2p1 p7p8 p2p1、 p7p8、p1p0 p8p0、p1p0 []
在构建该表之后,实际执行工作的代码只有几行。
注意:这里的代码使用复数值作为点。所以。real是。x。imag是。y。
def point_in_scantable(actives_table, events, xi, point):
beam = bisect(events, point.imag) - 1 # Binary search in sorted array.
actives_at_y = actives_table[beam]
total = sum([point.real > xi(e, point.imag) for e in actives_at_y])
return bool(total % 2)
我们对事件进行二进制搜索,以找到特定值的actives_at_y。对于在y点的所有活动,我们计算在我们点的特定y点的x段值。每次x截距大于点的x分量时加1。然后对总数乘以2。(这是偶数-奇数填充规则,你可以很容易地适应任何其他填充规则)。
完整的代码:
from bisect import bisect
def build_edge_list(polygon):
edge_list = []
for i in range(1, len(polygon)):
if (polygon[i].imag, polygon[i].real) < (polygon[i - 1].imag, polygon[i - 1].real):
edge_list.append((polygon[i], i))
edge_list.append((polygon[i - 1], ~i))
else:
edge_list.append((polygon[i], ~i))
edge_list.append((polygon[i - 1], i))
def sort_key(e):
return e[0].imag, e[0].real, ~e[1]
edge_list.sort(key=sort_key)
return edge_list
def build_scanbeam(edge_list):
actives_table = []
events = []
y = -float("inf")
actives = []
for pt, index in edge_list:
if y != pt.imag:
actives_table.append(list(actives))
events.append(y)
if index >= 0:
actives.append(index)
else:
actives.remove(~index)
y = pt.imag
return actives_table, events
def point_in_polygon(polygon, point):
def x_intercept(e, y):
pt0 = polygon[e-1]
pt1 = polygon[e]
if pt1.real - pt0.real == 0:
return pt0.real
m = (pt1.imag - pt0.imag) / (pt1.real - pt0.real)
b = pt0.imag - (m * pt0.real)
return (y - b) / m
edge_list = build_edge_list(polygon)
actives_table, events = build_scanbeam(edge_list)
try:
if len(point):
return [point_in_scantable(actives_table, events, x_intercept, p) for p in point]
except TypeError:
return point_in_scantable(actives_table, events, x_intercept, point)
def point_in_scantable(actives_table, events, xi, point):
beam = bisect(events, point.imag) - 1 # Binary search in sorted array.
actives_at_y = actives_table[beam]
total = sum([point.real > xi(e, point.imag) for e in actives_at_y])
return bool(total % 2)
如果忽略,则扫描表的构建时间为O(n*log(n))。我们实际上是在O(log(n) + k)时间内查到的。其中n是多边形中段数的大小,k是该多边形中典型的活动边数。其他的光线追踪方法实际上需要O(n)时间。每次我们检查一个点,它迭代整个多边形。所以即使有这个明显的次优实现,它也轻而易举地打败了其他所有的。
There's a few performance tricks that could be done, for example, we can lower the time complexity to O(log(n) + log(k)) time. To do this we would implement Bentley-Ottmann into the sweep line, and rather than processing the intersections as different events, we split the lines at the intersections. We then also sort the active edges by their x-intercepts. We then know that no intersections occur during a scanbeam and since we sorted our segments (taking care to order them correctly within the scanbeam even if they start at the same initial point (you need to look at the slopes, or just compare midpoints of the segments). We then have a sorted intersection-less actives lists scanbeam table which means we can binary search into active edge list as well. Though that sounds like a lot of work for a value of k which is going to be typically 2 or maybe 4.
此外,由于这基本上变成了一个查找表和一些x截距的最小计算,它更能用GPU完成。你不再需要在多边形上循环了。所以你可以用numpy这样的东西来大量计算这些点,这样你就可以一次做所有的计算,从而提高性能。
这只适用于凸形状,但是Minkowski Portal Refinement和GJK也是测试一个点是否在多边形中的很好的选择。您使用闵可夫斯基减法从多边形中减去点,然后运行这些算法来查看多边形是否包含原点。
另外,有趣的是,你可以用支持函数更隐式地描述你的形状,它以一个方向向量作为输入,并输出沿该向量的最远点。这可以让你描述任何凸形状..弯曲的,由多边形制成的,或混合的您还可以执行一些操作,将简单支持函数的结果组合起来,以生成更复杂的形状。
更多信息: http://xenocollide.snethen.com/mpr2d.html
此外,game programming gems 7讨论了如何在3d中做到这一点(:
nirg回答的Swift版本:
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
}
}
在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倍)。
Like David Segonds' answer suggests I use an approach of angle summation derived from my concave polygon drawing algorithm. It relies of adding up the approximate angles of subtriangles around the point to obtain a weight. A weight around 1.0 means the point is inside the triangle, a weight around 0.0 means outside, a weight around -1.0 is what happens when inside the polygon but in reverse order (like with one of the halves of a bowtie-shaped tetragon) and a weight of NAN if exactly on an edge. The reason it's not slow is that angles don't need to be estimated accurately at all. Holes can be handled by treating them as separate polygons and subtracting the weights.
typedef struct { double x, y; } xy_t;
xy_t sub_xy(xy_t a, xy_t b)
{
a.x -= b.x;
a.y -= b.y;
return a;
}
double calc_sharp_subtriangle_pixel_weight(xy_t p0, xy_t p1)
{
xy_t rot, r0, r1;
double weight;
// Rotate points (unnormalised)
rot = sub_xy(p1, p0);
r0.x = rot.x*p0.y - rot.y*p0.x;
r0.y = rot.x*p0.x + rot.y*p0.y;
r1.y = rot.x*p1.x + rot.y*p1.y;
// Calc weight
weight = subtriangle_angle_approx(r1.y, r0.x) - subtriangle_angle_approx(r0.y, r0.x);
return weight;
}
double calc_sharp_polygon_pixel_weight(xy_t p, xy_t *corner, int corner_count)
{
int i;
xy_t p0, p1;
double weight = 0.;
p0 = sub_xy(corner[corner_count-1], p);
for (i=0; i < corner_count; i++)
{
// Transform corner coordinates
p1 = sub_xy(corner[i], p);
// Calculate weight for each subtriangle
weight += calc_sharp_subtriangle_pixel_weight(p0, p1);
p0 = p1;
}
return weight;
}
因此,对于多边形的每一段,都形成一个子三角形,并计算点,然后旋转每个子三角形以计算其近似角度并添加到权重。
调用subtriangle_angle_approx(y, x)可以替换为atan2(y, x) / (2.*pi),但是一个非常粗略的近似值就足够精确了:
double subtriangle_angle_approx(double y, double x)
{
double angle, d;
int obtuse;
if (x == 0.)
return NAN;
obtuse = fabs(y) > fabs(x);
if (obtuse)
swap_double(&y, &x);
// Core of the approximation, a very loosely approximate atan(y/x) / (2.*pi) over ]-1 , 1[
d = y / x;
angle = 0.13185 * d;
if (obtuse)
angle = sign(d)*0.25 - angle;
return angle;
}
推荐文章
- 如何找到Java堆大小和内存使用(Linux)?
- 如何在HTML5中改变视频的播放速度?
- 我如何提高ASP。NET MVC应用程序性能?
- 列表推导式和函数式函数比for循环更快吗?
- 圆线段碰撞检测算法?
- Pandas loc vs iloc vs at vs iat?
- 当WebSockets可用时,为什么要使用AJAX ?
- SurfaceView和View的区别?
- 一个字符串的字符串列表
- .NET反射的成本有多高?
- 在c#中检查字符串是否只包含数字的最快方法
- JavaScript变量声明在循环外还是循环内?
- 为什么C语言这么快,为什么其他语言没有这么快或更快?
- String与StringBuilder
- 为什么引入无用的MOV指令会加速x86_64汇编中的紧循环?