- UID
- 5280
- 积分
- 9466
- 精华
- 贡献
-
- 威望
-
- 活跃度
-
- D豆
-
- 在线时间
- 小时
- 注册时间
- 2002-5-18
- 最后登录
- 1970-1-1
|
马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？立即注册
×
Determining the inclusion of a point P in a 2D planar polygon is a geometric problem that results in interesting algorithms. Two commonly used methods are: - The Crossing Number (cn) method
- which counts the number of times a ray starting from the point P crosses the polygon boundary edges. The point is outside when this "crossing number" is even; otherwise, when it is odd, the point is inside. This method is sometimes referred to as the "even-odd" test. - The Winding Number (wn) method
- which counts the number of times the polygon winds around the point P. The point is outside only when this "winding number" wn = 0; otherwise, the point is inside.
If a polygon is simple (i.e., it has no self intersections), then both methods give the same result for all points. But for non-simple polygons, the two methods can give different answers for some points. For example, when a polygon overlaps with itself, then points in the region of overlap are found to be outside using the crossing number, but are inside using the winding number, as shown in the diagrams: |
| | Crossing Number Method | | | |
| |
|
| | Vertices are numbered: 0 1 2 3 4 5 6 7 8 9 |
|
In this example, points inside the overlap region have wn = 2, implying that they are inside the polygon twice. Clearly, the winding number gives a better intuitive answer than the crossing number does. Despite this, the crossing number method is more commonly used since cn is erroneously thought to be significantly (up to 20 times!) more efficient to compute than wn [O'Rourke, 1998]. But this is not the case, and wn can be computed with the same efficiency as cn by counting signed crossings. In fact, our implementation for wn_PnPoly() is faster than the code for cn given by [Franklin, 2000], [ Haines, 1994] or [O'Rourke, 1998], although the cn "PointInPolygon()" routine of [Moller & Haines, 1999] is comparable to our wn algorithm. But the bottom line is that for both geometric correctness and efficiency reasons, the wn algorithm should always be preferred for determining the inclusion of a point in a polygon.
The Crossing NumberThis method counts the number of times a ray starting from a point P crosses a polygon boundary edge separating it's inside and outside. If this number is even, then the point is outside; otherwise, when the crossing number is odd, the point is inside. This is easy to understand intuitively. Each time the ray crosses a polygon edge, its in-out parity changes (since a boundary always separates inside from outside, right?). Eventually, any ray must end up beyond and outside the bounded polygon. So, if the point is inside, the sequence of crossings ">" must be: in > out > ... > in > out, and there are an odd number of them. Similarly, if the point is outside, there are an even number of crossings in the sequence: out > in > ... > in > out.
In implementing an algorithm for the cn method, one must insure that only crossings that change the in-out parity are counted. In particular, special cases where the ray passes through a vertex must be handled properly. These include the following types of ray crossings: Further, one must decide whether a point on the polygon's boundary is inside or outside. A standard convention is to say that a point on a left or bottom edge is inside, and a point on a right or top edge is outside. This way, if two distinct polygons share a common boundary segment, then a point on that segment will be in one polygon or the other, but not both at the same time. This avoids a number of problems that might occur, especially in computer graphics displays. A straightforward "crossing number" algorithm selects a horizontal ray extending to the right of P and parallel to the positive x-axis. Using this specific ray, it is easy to compute the intersection of a polygon edge with it. It is even easier to determine when no such intersection is possible. To count the total crossings, cn, the algorithm simply loops through all edges of the polygon, tests each for a crossing, and increments cn when one occurs. Additionally, the crossing tests must handle the special cases and points on an edge. This is accomplished by the
Edge Crossing Rules- an upward edge includes its starting endpoint, and excludes its final endpoint;
- a downward edge excludes its starting endpoint, and includes its final endpoint;
- horizontal edges are excluded
- the edge-ray intersection point must be strictly right of the point P.
One can apply these rules to the preceding special cases, and see that they correctly determine valid crossings. Note that Rule #4 results in points on a right-side boundary edge being outside, and ones on a left-side edge being inside.
Pseudo-Code: Crossing # InclusionCode for this algorithm is well-known, and the edge crossing rules are easily expressed. For a polygon represented as an array V[n+1] of vertex points with
, popular implementation logic ([Franklin, 2000], [O'Rourke, 1998]) is as follows: [C++] 纯文本查看 复制代码 typedef struct {int x, y;} Point;
cn_PnPoly( Point P, Point V[], int n )
{
int cn = 0; // the crossing number counter
// loop through all edges of the polygon
for (each edge E[i]:V[i]V[i+1] of the polygon) {
if (E[i] crosses upward ala Rule #1
|| E[i] crosses downward ala Rule #2) {
if (P.x < x_intersect of E[i] with y=P.y) // Rule #4
++cn; // a valid crossing to the right of P.x
}
}
return (cn&1); // 0 if even (out), and 1 if odd (in)
} |
|
Note that the tests for upward and downward crossing satisfying Rules #1 and #2 also exclude horizontal edges (Rule #3). All-in-all, a lot of work is done by just a few tests which makes this an elegant algorithm. However, the validity of the crossing number method is based on the "Jordan Curve Theorem" which says that a simple closed curve divides the 2D plane into exactly 2 connected components: a bounded "inside" one and an unbounded "outside" one. The catch is that the curve must be simple (without self intersections), otherwise there can be more than 2 components and then there is no guarantee that crossing a boundary changes the in-out parity. For a closed (and thus bounded) curve, there is exactly one unbounded "outside" component; but bounded components can be either inside or outside. And two of the bounded components that have a shared boundary may both be inside, and crossing over their shared boundary would not change the in-out parity.
The Winding NumberOn the other hand, the winding number accurately determines if a point is inside a nonsimple closed polygon. It does this by computing how many times the polygon winds around the point. A point is outside only when the polygon doesn't wind around the point at all which is when the winding number wn = 0. More generally, one can define the winding number
of any closed continuous curve C around a point P in the 2D plane. Let the continuous 2D curve C be defined by the points
, for
with
. And let P be a point not on C. Then, define the vector c(P,u) = C(u) – P from P to C(u), and the unit vector w(P,u) = c(P,u) / |c(P,u)| which gives a continuous function
mapping the point C(u) on C to the point w(P,u) on the unit circle
. This map can be represented in polar coordinates as
where is a positive counterclockwise angle in radians. The winding number is then equal to the integer number of times that W(P) wraps C around S1. This corresponds to a homotopy class of S1, and can be computed by the integral:
When the curve C is a polygon with vertices V0,V1,...,Vn = V0, this integral reduces to the sum of the (signed) angles that each edge ViVi+1 subtends with the point P. So, if
, we have:
This formula is clearly not very efficient since it uses a computationally expensive arccos() trig function. But, a simple observation lets us replace this formula by a more efficient one. Pick any point Q on S1. Then, as the curve W(P) wraps around S1, it passes Q a certain number of times. If we count (+1) when it passes Q counterclockwise, and (–1) when it passes clockwise, then the accumulated sum is exactly the total number of times that W(P) wraps around S1, and is equal to the winding number . Further, if we take an infinite ray R starting at P and extending in the direction of the vector Q, then intersections where R crosses the curve C correspond to the points where W(P) passes Q. To do the math, we have to distinguish between positive and negative crossings where C crosses R from right-to-left or left-to-right. This can be determined by the sign of the dot product between a normal vector to C and the direction vector q = Q [Foley et al, 1996, p-965], and when the curve C is a polygon, one just needs to make this determination once for each edge. For a horizontal ray R from P, testing whether an edge's endpoints are above and below the ray suffices. If the edge crosses the positive ray from below to above, the crossing is positive (+1); but if it crosses from above to below, the crossing is negative (–1). One then simply adds all crossing values to get . For example:
Additionally, one can avoid computing the actual edge-ray intersection point by using the isLeft() attribute; however, it needs to be applied differently for ascending and descending edges. If an upward edge crosses the ray to the right of P, then P is on the left side of the edge since the triangle ViVi+1P is oriented counterclockwise. On the other hand, a downward edge crossing the positive ray would have P on the right side since the triangle ViVi+1P would then be oriented clockwise.
Pseudo-Code: Winding Number InclusionThis results in the following wn algorithm which is an adaptation of the cn algorithm and uses the same edge crossing rules as before to handle special cases. [C++] 纯文本查看 复制代码 typedef struct {int x, y;} Point;
wn_PnPoly( Point P, Point V[], int n )
{
int wn = 0; // the winding number counter
// loop through all edges of the polygon
for (each edge E[i]:V[i]V[i+1] of the polygon) {
if (E[i] crosses upward ala Rule #1) {
if (P is strictly left of E[i]) // Rule #4
++wn; // a valid up intersect right of P.x
}
else
if (E[i] crosses downward ala Rule #2) {
if (P is strictly right of E[i]) // Rule #4
--wn; // a valid down intersect right of P.x
}
}
return wn; // =0 <=> P is outside the polygon
} |
|
Clearly, this winding number algorithm has the same efficiency as the analogous crossing number algorithm. Thus, since it is more accurate in general, the winding number algorithm should always be the preferred method to determine inclusion of a point in an arbitrary polygon. The wn algorithm's efficiency can be improved further by rearranging the crossing comparison tests. This is shown in the detailed implementation of wn_PnPoly() given below. In that code, all edges that are totally above or totally below P get rejected after only two (2) inequality tests. However, currently popular implementations of the cn algorithm ([Franklin, 2000], [ Haines, 1994], [O'Rourke, 1998]) use at least three (3) inequality tests for each rejected edge. Since most of the edges in a large polygon get rejected in practical applications, there is about a 33% (or more) reduction in the number of comparisons done. In runtime tests using very large (1,000,000 edge) random polygons (with edge length < 1/10 the polygon diameter) and 1000 random test points (inside the polygon's bounding box), we measured a 20% increase in efficiency overall.
EnhancementsThere are some enhancements to point in polygon algorithms [Haines, 1994] that software developers should be aware of. We mention a few that pertain to ray crossing algorithms. However, there are other techniques that give better performance in special cases such as testing inclusion in small convex polygons like triangles. These are discussed in [Haines, 1994].
Bounding Box or BallIt is efficient to first test that a point P is inside the bounding box or ball of a large polygon before testing all edges for ray crossings. If a point is outside the bounding box or ball, it is also outside the polygon, and no further testing is needed. But, one must precompute the bounding box (the max and min for vertex x and y coordinates) or the bounding ball (center and minimum radius) and store it for future use. This is worth doing if more than a few points are going to be tested for inclusion, which is generally the case. Further information about computing bounding containers can be found in Algorithm 8.
3D Planar PolygonsIn 3D applications, one sometimes wants to test a point and polygon that are in the same plane. For example, one may have the intersection point of a ray with the plane of a polyhedron's face, and want to test if it is inside the face. Or one may want to know if the base of a 3D perpendicular dropped from a point is inside a planar polygon. 3D inclusion is easily determined by projecting the point and polygon into 2D. To do this, one simply ignores one of the 3D coordinates and uses the other two. To optimally select the coordinate to ignore, compute a normal vector to the plane, and select the coordinate with the largest absolute value [Snyder & Barr, 1987]. This gives the projection of the polygon with maximum area, and results in robust computations.
Convex or Monotone PolygonsWhen a polygon is known to be convex, many computations can be speeded up. For example, the bounding box can be computed in O(log n) time [O'Rourke, 1998] instead of O(n) as for nonconvex polygons. Also, point inclusion can be tested in O(log n) time. More generally, the following algorithm works for convex or monotone polygons. A convex or y-monotone polygon can be split into two polylines, one with edges increasing in y, and one with edges decreasing in y. The split occurs at two vertices with max y and min y coordinates. Note that these vertices are at the top and bottom of the polygon's bounding box, and if they are both above or both below P’s y-coordinate, then the test point is outside the polygon. Otherwise, each of these two polylines intersects a ray paralell to the x-axis once, and each potential crossing edge can be found by doing a binary search on the vertices of each polyline. This results in a practical O(log n) (preprocessing and runtime) algorithm for convex and montone polygon inclusion testing. For example, in the following diagram, only three binary search vertices have to be tested on each polyline (A1, A2, A3 ascending; and B1, B2, B3 descending): This method also works for polygons that are monotone in an arbitrary direction. One then uses a ray from P that is perpendicular to the direction of monotonicity, and the algorithm is easily adapted. A little more computation is needed to determine which side of the ray a vertex is on, but for a large enough polygon, the O(log n) performance more than makes up for the overhead.
ImplementationsHere is a "C++" implementation of the winding number algorithm for the inclusion of a point in polygon. We just give the 2D case, and use the simplest structures for a point and a polygon which may differ in your application. [C++] 纯文本查看 复制代码 // Copyright 2000 softSurfer, 2012 Dan Sunday
// This code may be freely used and modified for any purpose
// providing that this copyright notice is included with it.
// SoftSurfer makes no warranty for this code, and cannot be held
// liable for any real or imagined damage resulting from its use.
// Users of this code must verify correctness for their application.
// a Point is defined by its coordinates {int x, y;}
//===================================================================
// isLeft(): tests if a point is Left|On|Right of an infinite line.
// Input: three points P0, P1, and P2
// Return: >0 for P2 left of the line through P0 and P1
// =0 for P2 on the line
// <0 for P2 right of the line
// See: Algorithm 1 "Area of Triangles and Polygons"
inline int
isLeft( Point P0, Point P1, Point P2 )
{
return ( (P1.x - P0.x) * (P2.y - P0.y)
- (P2.x - P0.x) * (P1.y - P0.y) );
}
//===================================================================
// cn_PnPoly(): crossing number test for a point in a polygon
// Input: P = a point,
// V[] = vertex points of a polygon V[n+1] with V[n]=V[0]
// Return: 0 = outside, 1 = inside
// This code is patterned after [Franklin, 2000]
int
cn_PnPoly( Point P, Point* V, int n )
{
int cn = 0; // the crossing number counter
// loop through all edges of the polygon
for (int i=0; i<n; i++) { // edge from V[i] to V[i+1]
if (((V[i].y <= P.y) && (V[i+1].y > P.y)) // an upward crossing
|| ((V[i].y > P.y) && (V[i+1].y <= P.y))) { // a downward crossing
// compute the actual edge-ray intersect x-coordinate
float vt = (float)(P.y - V[i].y) / (V[i+1].y - V[i].y);
if (P.x < V[i].x + vt * (V[i+1].x - V[i].x)) // P.x < intersect
++cn; // a valid crossing of y=P.y right of P.x
}
}
return (cn&1); // 0 if even (out), and 1 if odd (in)
}
//===================================================================
// wn_PnPoly(): winding number test for a point in a polygon
// Input: P = a point,
// V[] = vertex points of a polygon V[n+1] with V[n]=V[0]
// Return: wn = the winding number (=0 only when P is outside)
int
wn_PnPoly( Point P, Point* V, int n )
{
int wn = 0; // the winding number counter
// loop through all edges of the polygon
for (int i=0; i<n; i++) { // edge from V[i] to V[i+1]
if (V[i].y <= P.y) { // start y <= P.y
if (V[i+1].y > P.y) // an upward crossing
if (isLeft( V[i], V[i+1], P) > 0) // P left of edge
++wn; // have a valid up intersect
}
else { // start y > P.y (no test needed)
if (V[i+1].y <= P.y) // a downward crossing
if (isLeft( V[i], V[i+1], P) < 0) // P right of edge
--wn; // have a valid down intersect
}
}
return wn;
}
//===================================================================
下面是LISP实现代码：
;;=========================================================================; ;; ; ;; PtInPoly_p(): Winding number test for a point in a polygon ; ;; ; ;; Input: p = a point ; ;; en = ename or vla-object of a closed polyline. ; ;; ; ;; Return: t, if point is in polyline. (wn is not 0) ; ;; ; ;; Original code in C++ by Dan Sunday, Translated by ymg ; ;;=========================================================================; (defun ptinpoly_p (p l / i p1 p2 v wn y isleft) ;; Returns t if point p is strictly to the left of v1->v2 by ymg ; (defun isleft (p v1 v2 / xp yp) (setq xp (car p) yp (cadr p) ) (- (* (- (cadr v1) yp) (- (car v2) xp)) (* (- (car v1) xp) (- (cadr v2) yp))) ) ;; Initialize the loop. ; (setq p1 (car l) y (cadr p) wn 0 ) ;; Loop through all edges of polygon. ; (foreach p2 (cdr l) (if (<= (cadr p1) y) (if (> (cadr p2) y) (if (> (isleft p p1 p2) 0) (setq wn (1+ wn)) ) ) (if (<= (cadr p2) y) (if (< (isleft p p1 p2) 0) (setq wn (1- wn)) ) ) ) (setq p1 p2) ) (not (zerop wn)) ;; Returns t if wn not equal to zero. ; ) (princ "\nSelect Points: ") (if (setq ss (ssget '((0 . "point") (8 . "0,testpoint")))) (progn (foreach e (vl-remove-if 'listp (mapcar 'cadr (ssnamex ss))) (if (= "0" (cdr (assoc 8 (entget e)))) (progn (setq b1 (cons (setq p (cdr (assoc 10 (entget e)))) b1)) (entmakex (list '(0 . "POINT") '(8 . "0") (cons 10 (list (+ 100 (car p)) (cadr p) (last p)))) ) ) (progn (setq tp1 (cons (cdr (assoc 10 (entget e))) tp1)) (entdel e)) ) ) (setq b2 (mapcar '(lambda (p) (list (+ 100 (car p)) (cadr p) (last p))) b1)) (setq tp2 (mapcar '(lambda (p) (list (+ 100 (car p)) (cadr p) (last p))) tp1)) (foreach p tp1 (if (ptinpoly_p p b1) (entmakex (list '(0 . "POINT") (cons 10 p) '(62 . 1))) (entmakex (list '(0 . "POINT") (cons 10 p) '(62 . 7))) ) ) (foreach p tp2 (if (inside_p p b2) (entmakex (list '(0 . "POINT") (cons 10 p) '(62 . 1))) (entmakex (list '(0 . "POINT") (cons 10 p) '(62 . 7))) ) ) ) ) (princ) ) Tomas Moller & Eric Haines, "Ray/Polygon Intersection" in Real-Time Rendering (3rd Edition) (2008) John M. Snyder & Alan H. Barr, "Ray Tracing Complex Models Containing Surface Tessellations", Computer Graphics 21(4), 119-126 (1987) [also in the Proceedings of SIGGRAPH 1987] |
-
|
如果你在论坛求助问题，并且已经从坛友或者管理的回复中解决了问题，请把帖子标题加上【已解决】；
如何回报帮助你解决问题的坛友，一个好办法就是给对方加【D豆】，加分不会扣除自己的积分，做一个热心并受欢迎的人！
|申请友链|Archiver|手机版|小黑屋|辽公网安备|晓东CAD家园
( 辽ICP备15016793号 )
窥视卡
雷达卡
已领礼包: 
财富等级: 
发表于 2017-7-14 17:57:46
提升卡
置顶卡
沉默卡
变色卡
千斤顶
发表于 2017-7-15 15:35:14
