多边形封闭区域算法

前些天在考虑一个几何算法,关于如何判断一个点是否在一个给定的多边形内部。这应该是一个比较常规的算法,我以前对几何算法了解的不多,所以既然想到了就稍微研究了一下。
查了一下相关的资料,目前有几个O(N)的算法,其中N是多边形的顶点数。

第一个叫做交替(Alternative)算法。 如下图所示

交替法


算法检查所给定的点和在该点右边的多边形区域的边界的交点数,若交点数为奇数,则该点在多边形内部,若交点数为偶数,则该点在多边形外部。可以想象成从一个点向右发出一条射线,然后计算这条射线与多边形边界的交点数。从上图中标出的蓝色点向右射出的射线与多边形有1个交点,故判定其在多边形内部,而从标出的白色点向右射出的射线则有偶数个交点,故判定其在多边形外部。

这个算法的叙述起来很平常,但实现起来却有两个略带Tricky的地方。

其一是坐标精度引起的,因为计算机对于坐标的描述是离散的,在计算向右射出的射线与多边形交点的时候很容易出现射线与多边形的某几条边平行的情况,对于这种情况,我们的算法需要向前多看几个点,以确定射线是确实与多边形相交还是只是“擦过”而已。

其二是对于非常靠近多边形边界的点的判定,你可以认为这些点是在多边形外,也可以认为是在多边形内,但你必须在算法中描述这些情况。

感兴趣的可以了解一下我的算法实现,附在文末。

这个算法在多边形为简单多边形的情况下是没有异议的,但在复杂多边形的情况下会有一些疑问,这涉及到“在多边形内还是多边形外”的定义问题。如下图

2009-09-24_122701

图中,如果用Alternative算法,则被多边形包围的那个点被判断为在多边形外部。

如果你认为这种定义不符合你的需要,那么有另外一种算法,称为Winding算法

Winding算法的思路是这样的,将你置身于所给定的点上,然后让你的视线从多边形边界上的一点开始,选择一个方向绕着多边形的边界转一圈,如果你发现你的身体在这个过程中转了360度的非零整数倍,那么你所站在的这个点在多边形的内部,如果你的身体并没有转动(事实上你转动了,只是正向和逆向的转动抵消了),那么就在多边形的外部。

这个算法我没有去实现,一般来说它要比Alternative算法慢,因为它需要计算每条多边形的边与多给定的点之间的夹角。

关于算法实现:

程序的界面是用WTL写的,如果你想要编译源代码,请在Visual Studio中指定WTL的include文件夹位置。

算法核心代码如下,其中vSelectionPoints是多边形

C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
// thePoint:给定的点
// vSelectionPoints:给定的多边形顶点向量
BOOL CheckInside(CPoint thePoint, std::vector vSelectionPoints)
{
    int nCrossCount = 0;
    for (size_t j = 0; j < vSelectionPoints.size(); ++j)
    {
        CPoint &p1 = vSelectionPoints[j];
        CPoint &p2 = vSelectionPoints[(j+1)%vSelectionPoints.size()];
        if (thePoint == p1 || thePoint == p2)
        {
            nCrossCount = 0;
            break;
        }
        else if (thePoint.y > max(p1.y, p2.y) || thePoint.y < min(p1.y, p2.y))
        {
            continue;
        }
        else if (thePoint.y < max(p1.y, p2.y) && thePoint.y > min(p1.y, p2.y))
        {
            if (p2.x == p1.x)
            {
                if (thePoint.x == p1.x)
                {
                    nCrossCount = 0;
                    break;
                }
                else if (thePoint.x < p1.x)
                {
                    nCrossCount++;
                }
            }
            else
            {
                double fCrossX = (thePoint.y - p1.y) * double(p2.x - p1.x) / double(p2.y - p1.y) + p1.x;
                if (fabs(fCrossX - thePoint.x) < ESPILON)
                {
                    nCrossCount = 0;
                    break;
                }
                if (fCrossX > thePoint.x)
                {
                    nCrossCount++;
                }
            }
        }
        else
        {
            if (thePoint.y != p1.y && thePoint.y == p2.y && thePoint.x < p2.x)
            {
                while(TRUE)
                {
                    j++;
                    CPoint *pp3 = &vSelectionPoints[(j+1)%vSelectionPoints.size()];
                    if (pp3->y != p2.y)
                    {
                        break;
                    }
                }
 
                CPoint &p3 = vSelectionPoints[(j+1)%vSelectionPoints.size()];
 
                if (thePoint.y <= max(p1.y, p3.y) && thePoint.y >= min(p1.y, p3.y))
                {
                    nCrossCount++;
                }
            }
        }
    }
    if (nCrossCount % 2 != 0)
    {
        return TRUE;
    }
    return FALSE;
}

工程文件下载:
WTLSelectionWindow

Comments (16)

  1. 10:43, 2009-09-25boolzhu  / Reply

    恩,我写过这个解法的三维版本,判断一个点是否在多边形内部,但问题很多,现在只是勉强能用。最大问题就是奇偶性常常会在面片的相交处出错,会在重合的点或者线上多算相交的次数。你不如写个三维版本给我用吧,哈哈~

  2. 14:19, 2009-09-25卢松松  / Reply

    你还好学术研究了?

  3. 14:42, 2009-09-25Justice  / Reply

    @卢松松
    我们这里是3个人的额...
    研究算法的是大牛Chris童鞋 = =#

  4. 19:46, 2009-09-25shen  / Reply

    大顶技术贴!!!经小鹤鉴定,该算法很有用,值得研究~~期待有三维的版本~~~

  5. 20:23, 2009-09-25Chris  / Reply

    @boolzhu
    3D的话不知道怎么做成GUI= =b

  6. 23:12, 2009-09-25Sinnyn  / Reply

    @Chris
    你做个实现给他个interface

  7. 20:40, 2009-09-26StepInto  / Reply

    赞WTL!!
    这个Winding算法有个改进版,就是用cross product求带符号面积,Sigma area(Ai, P, Ai+1)
    如查面积和是0(正负相抵)那么在形外,否则结果应该是多边形的面积就是在形内

  8. 22:34, 2009-09-26Chris  / Reply

    @StepInto
    恩,好办法,这样就可以省去求arccos或者arctan的计算,而且可以用在3D的情况下

  9. 23:55, 2009-09-26Melvin Xie  / Reply

    让我想起了何媛军的图形学~赞

  10. 10:59, 2009-09-27Chris  / Reply

    @Melvin Xie
    啊,惊现谢mm~

  11. 10:07, 2009-10-24hayate  / Reply

    好奇的问,这个demo为啥要用华丽丽的WTL呢,另外窗口绘图是用GDI么,还是别的什么

  12. 10:50, 2009-10-24hayate  / Reply

    @hayate
    好吧,我已经看过source乐。。。
    感谢,很有趣的程序

  13. 13:49, 2009-10-24Chris  / Reply

    @hayate
    其实我是在学WTL= =,想到了这个算法,就顺便实现一下,哈

  14. 00:15, 2009-10-26半瓶墨水  / Reply

    图文并茂的描述很赞

    但是代码加色在RSS阅读器里面表现很差 - 白色的字体配上白背景

    呵呵,欢迎试用代码发芽网的语法高亮 http://fayaa.com/code/

  15. 09:59, 2009-10-26Chris  / Reply

    @半瓶墨水
    广告贴~呵呵,好,我试试看

  16. 15:04, 2009-10-26Justice  / Reply

    @半瓶墨水
    你们fayaa的我试过,但是调来调去感觉效果还是WP-Syntax弄出来舒服呀,而且这个是自动的~你说的那个问题也可以解决的,加个inline的background就可以了~

Leave a Reply

Allowed Tags - You may use these HTML tags and attributes in your comment.

<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong> <pre lang="" line="" escaped="" highlight="">

Pingbacks (0)

› No pingbacks yet.