几种判断点与多边形关系的算法介绍

本文讨论如何判断一个点是在多边形内部，边上还是在外部。为了方便，这里的多边形默认为有向多边形，规定沿多边形的正向，边的左侧为多边形的内侧域，即多边形边按逆时针方向遍历，不考虑自交等复杂情况。

比较常见的判断点与多边形关系的算法有射线法、面积法、点线判断法和弧长法等，算法复杂度都为O(n)，不过只有射线法可以正确用于凹多边形，其他3个只可以用于凸多边形。

1. 射线法

射线法是使用最广泛的算法，这是由于相比较其他算法而言，它不但可以正确使用在凹多边形上，而且不需要考虑精度误差问题。该算法思想是从点出发向右水平做一条射线，计算该射线与多边形的边的相交点个数，当点不在多边形边上时，如果是奇数，那么点就一定在多边形内部，否则，在外部。

typedef struct{
    double x;
    double y;
}point;
typedef struct{
    int numpoint;
    point *pointlist; 
}polygon;
int raycasting(polygon *pol, point *ref){
    int count = 0;
    for(int i=0; i<pol->numpoint; ++i){
        point *p1 = pol->pointlist+i;
        point *p2 = i==pol->numpoint-1?pol->pointlist:pol->pointlist+i+1;
        if((ref->y>=p1->y&&ref->y<=p2->y) || (ref->y>=p2->y&&ref->y<=p1->y)){
            double t = (ref->y-p1->y)/(p2->y-p1->y);
            double xt = p1->x+t*(p2->x-p1->x);
            if(ref->x==xt) return ONSIDE;
            if(ref->x<xt) ++count;
        }
    }
    return count%2?INSIDE:OUTSIDE;
}

2. 面积法

面积法的思想是如果点在多边形内部或者边上，那么点与多边形所有边组成的三角形面积和等于多边形面积。多边形的面积公式可以用叉积计算，在上一篇已经介绍了。不过计算面积是会有一定误差的，需要设置精度的误差范围。

int areacheck(polygon *pol, point *ref){
    double polygonarea = calpolygonarea(pol);
    double triarea = 0;
    for(int i=0; i<pol->numpoint; i++){
        point *p1 = pol->pointlist+i;
        point *p2 = i==pol->numpoint-1?pol->pointlist:pol->pointlist+i+1;
        double area = ABS(cross(ref, p1, p2));
        if(area==0) return ONSIDE;
        triarea += 0.5*area;
    }
    return ABS(polygonarea-triarea)<TOL?INSIDE:OUTSIDE;
}

3. 点线判断法

对于多边形，如果一个点它的所有边的左边，那么这个点一定在多边形内部。利用叉积正好可以判断点与给定边的关系，即点是在边的左边右边还是边上。

int crosscheck(polygon *pol, point *ref){
    for(int i=0; i<pol->numpoint; i++){
        point *p1 = pol->pointlist+i;
        point *p2 = i==pol->numpoint-1?pol->pointlist:pol->pointlist+i+1;
        double cr = cross(ref, p1, p2);
        if(cr>0) return OUTSIDE;
    if(cr==0) return ONSIDE;
    }
    return INSIDE;
}

4. 弧长法

弧长法是改进后的转角法，解决了传统转角法的精度问题。算法思想是，以被测点O为坐标原点，将平面划分为4个象限，对每个多边形顶点P[i]，计算其所在的象限，然后顺序访问多边形的各个顶点P[i]，分析P[i]和P[i+1]，有下列三种情况：

（1）P[i+1]在P[i]的下一象限。此时弧长和加π/2；

（2）P[i+1]在P[i]的上一象限。此时弧长和减π/2；

（3）P[i+1]在Pi的相对象限。利用叉积f=y[i+1]*x[i]-x[i+1]*y[i]计算Op[i]与Op[i+1]的关系，若f=0，Op[i]与Op[i+1]共线，点在多边形边上；若f<0，Op[i+1]在Op[i]逆时针方向，弧长和减π；若f>0，Op[i+1]在Op[i]顺时针方向，弧长和加π。

由于顶点在原点和坐标轴上时需要特殊处理，为了准确性，应该在每次计算顶点时都用叉积判断P是否在当前边上，另外将π用整数代替可以避免浮点数的精度误差。

int arclengthcheck(polygon *pol, point *ref){
    int arc = 0, q1, q2; 
    double x1, y1, x2, y2;
    for(int i=0; i<pol->numpoint; i++){
        point *p1 = pol->pointlist+i;
        point *p2 = i==pol->numpoint-1?pol->pointlist:pol->pointlist+i+1;
        double cr = cross(p1, ref, p2);
        if(cr==0) return ONSIDE; 
        x1 = p1->x-ref->x;
        y1 = p1->y-ref->y;
        x2 = p2->x-ref->x;
        y2 = p2->y-ref->y;
        q1 = x1>0?(y1>0?0:3):(y1>0?1:2);
        q2 = x2>0?(y2>0?0:3):(y2>0?1:2);
        int g = (q2-q1+4)%4;
        if(g==1) arc++;
        else if(g==3) arc--;
        else if(g==2) cr>0?(arc+=2):(arc-=2);  
    }
    return arc==0?OUTSIDE:INSIDE;
}

code:

https://github.com/coderkian/algorithm/blob/master/inpolygon.c