幾種判斷點與多邊形關系的算法介紹

本文討論如何判斷一個點是在多邊形內部，邊上還是在外部。為了方便，這里的多邊形默認為有向多邊形，規定沿多邊形的正向，邊的左側為多邊形的內側域，即多邊形邊按逆時針方向遍歷，不考慮自交等復雜情況。

比較常見的判斷點與多邊形關系的算法有射線法、面積法、點線判斷法和弧長法等，算法復雜度都為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