判断点在多边形内还是外部

判断点在多边形内部还是外部是一个基础但很重要的算法。
下面介绍两种方法：
一：射线法：
从点P向左作水平射线，如果P在多边形内部，那么这条射线与多边形的交点必为奇数；
                                           如果P在多边形外部，则交点个数必为偶数（0 included）。
顺序考虑多边形的每条边，便可求出交点总个数，O(n)。
对于边PQ
特殊情况：（1）如果射线穿过顶点P或Q，则会被重复计数。
                    （2）如果PQ为水平线段，则有可能与射线“重合”。
处理办法：
                     对于（1）如果穿过的顶点纵坐标是PQ中较小的则忽略
                     对于（2）直接忽略  ？ P点要么在PQ上，要么不在。只需在开始判断一下是否在PQ上。

具体实现代码：

复制代码

1. #include <stdio.h>
2. #include <math.h>
3. #include <iostream>
4. using namespace std;
5. #define inf 1e9
6. #define exp 1e-9
8. struct node
9. {
10.     double x;
11.     double y;
12. };
13. node nodes[105], p;
14. struct segment
15. {
16.     node st;
17.     node en;
18. };
19. int n;
20. //叉积 if > 0 p1在p2的顺时针方向
21. double multi(node p1, node p2, node p0)
22. {
23.     return (p1.x - p0.x) * (p2.y - p0.y) - (p1.y - p0.y) * (p2.x - p0.x);
24. }
25. //if叉积为零说明共线，如果该点在以side为对角线的矩形中说明p在side上
26. int isonline(node p, segment side)
27. {
28.     return (fabs(multi(p, side.st, side.en)) < exp &&
29.         min(side.st.x, side.en.x) <= p.x && p.x <= max(side.st.x, side.en.x) &&
30.         p.y >= min(side.st.y, side.en.y) && p.y <= max(side.st.y,side.en.y));
31. }
32. //线段相交，跨立试验（叉积判断），快速排斥试验
33. int intersect(segment l, segment r)
34. {
35.     return (   max(l.st.x, l.en.x) >= min(r.st.x, r.en.x)//快速排斥试验用来确定以两线段为对角线的矩形是否相交
36.             && max(r.st.x, r.en.x) >= min(l.st.x, l.en.x)
37.             && max(l.st.y, l.en.y) >= min(r.st.y, r.en.y)
38.             && max(r.st.y, r.en.y) >= min(l.st.y, l.en.y)
39.             && multi(r.st, l.en, l.st) * multi(r.en, l.en, l.st) <= 0
40.             && multi(l.st, r.en, r.st) * multi(l.en, r.en, r.st) <= 0);
41. }
42. int isinpolygon(node p)
43. {
44.     int i, count = 0;
45.     segment line, side;
46.     if(n == 1) return fabs(p.x - nodes[0].x) < exp && fabs(p.y - nodes[0].y) < exp;
47.     if(n == 2)
48.     {
49.         side.st = nodes[0];
50.         side.en = nodes[1];
51.         return isonline(p, side);
52.     }
53.     line.st = p;
54.     line.en.y = p.y;
55.     line.en.x = -inf;
56.     for( i = 0; i < n; i ++ )
57.     {
58.         side.st = nodes[i];
59.         side.en = nodes[(i + 1) % n];
60.         if(isonline(p,side)) return 1;// 如果p在side上
61.         if(fabs(side.st.y - side.en.y) < exp) continue;//如果是水平的则不考虑
62.         //两端点相交line，取y坐标大的以避免重复计数
63.         if(isonline(side.st, line) && side.st.y > side.en.y) count ++;
64.         else if(isonline(side.en, line) && side.en.y > side.st.y) count ++;
65.         else if(intersect(line, side)) count ++;
66.     }
67.     return count % 2;
68. }
69. int main()
70. {
71.     while(scanf("%d",&n) && n)
72.     {
73.        for(int i = 0; i < n; i ++)
74.            scanf("%lf%lf",&nodes[i].x, &nodes[i].y);
75.        scanf("%lf%lf",&p.x,&p.y);
76.        if(isinpolygon(p)) printf("Within\n");
77.        else printf("Outside\n");
78.     }
79.     return 0;
80. }

二：改进的弧长法：摘自网上
有关"弧长法"的介绍:"弧长法要求多边形是有向多边形,一般规定沿多边形的正向,边的左侧为多边形的内侧域.以被测点为圆心作单位圆,将全部有向边向单位圆作径向投影,并计算其中单位圆上弧长的代数和,若代数和为0,则点在多边形外部;若代数和为2π,则点在多边形内部;若代数和为π,则点在多边形上."
根据上面的介绍,其实弧长法就是转角法,但它的改进方法比较比较厉害:将坐标原点平移到被测点P,这个新坐标系将平面划分为4个象限,对每个多边形顶点P,只考虑其所在的象限，然后按邻接顺序访问多边形的各个顶点P,分析P和P[i+1],有下列三种情况:
 (1) P[i+1]在P的下一象限,此时弧长和加π/2;
（2）P[i+1]在P的上一象限,此时弧长和减π/2;
（3）P[i+1]在P的相对象限,首先计算f=p[i+1].y*p.x-p[i+1].x*p.y(叉积),若f=0,则点在多边形上;若f<0,弧长和减π;若f>0,弧长和加π.最后对算出的代数和和上述的情况一样判断即可.
实现的时候还有两点要注意:
1> 若P的某个坐标为0时,一律当正号处理;
 2> 若被测点和多边形的顶点重合时要特殊处理.
还有一个问题那就是当多边形的某条边在坐标轴上而且两个顶点分别在原点的两侧时会出错,如边(3,0)……>(-3,0),按以上的处理,象限分别是第一和第二,这样会使代数和加π/2,有可能导致最后结果是被测点在多边形外,而实际上被测点是在多边形上(该边穿过该点).
对于这点,处理办法是:每次计算P和P[i+1]时,就计算叉积和点积,判断该点是否在该边上,是则判断结束,否则继续上述过程,这样牺牲了时间,但保证了正确性.
具体实现的时候,由于只需知道当前点和上一点的象限位置,所以附加空间只需O(1).实现的时候可以把上述的"π/2"改成1,"π"改成2,这样便可以完全使用整数进行计算,不必考虑顶点的顺序,逆时针和顺时针都可以处理,只是最后的代数和符号不同而已.

优点：实现简单，具有很高的精度，只做乘法和减法。

具体实现代码：

复制代码

1. #include <stdio.h>
2. #include <iostream>
3. using namespace std;
4. #define M 105
5. struct node
6. {
7.     double x, y;
8. }p[M];
10. int inpolygon(node t, int n)
11. {
12.     int i, t1, t2, sum, f;
13.     for(i = 0; i <= n; i ++) p[i].x -= t.x, p[i].y -= t.y;
14.     t1 = p[0].x>=0 ?(p[0].y>=0?0:3) :(p[0].y>=0?1:2);  
15.     for(sum = 0, i = 1; i <= n; i ++)
16.     {
17.         if(!p[i].x && !p[i].y) break;            
18.         f = p[i].y * p[i-1].x - p[i].x * p[i-1].y;  
19.         if(!f && p[i-1].x*p[i].x <= 0 && p[i-1].y*p[i].y <= 0) break;  
20.         t2 = p[i].x>=0 ?(p[i].y>=0?0:3) :(p[i].y>=0?1:2);  
21.         if(t2 ==(t1 + 1) % 4) sum += 1;        
22.         else if(t2 ==(t1 + 3) % 4) sum -= 1;  
23.         else if(t2 ==(t1 + 2) % 4)              
24.         {                                            
25.             if(f > 0) sum += 2;
26.             else sum -= 2;
27.         }
28.         t1 = t2;
29.     }
30.     if(i<=n || sum) return 1;
31.     return 0;
32. }
33. int main()
34. {
35.     int n, i;
36.     node t;
37.     while(scanf("%d", &n) && n)
38.     {
39.         for(i = 0; i < n; i ++)
40.             scanf("%lf%lf", &p[i].x, &p[i].y);
41.         p[n] = p[0];
42.         scanf("%lf%lf", &t.x, &t.y);
43.         if(inpolygon(t, n)) printf("inside\n");
44.         else printf("outside\n");
45.     }
46.     return 0;
47. }

[ 此帖被高尚博在2011-08-12 11:22重新编辑 ]