判斷點在多邊形內部

0.前言

最近不斷遇到類似的幾何位置問題，一直沒有花時間去總結，本文總結了我常用點跟多邊形的位置判斷方法以及代碼。希望能夠對大家有所協助。

文中所指的多邊形均為凸多邊形，一些描述可能有誤，歡迎指正。

1.測試的多邊形

在開始之前，我們需要先構建好測試環境。

我構建了一個比較特殊的多邊形，如下。

/ \

|  |

|_|

從最上面的頂點順時針座標（螢幕座標系）分別為：(40,10) (60,30)(60,50) (20,50) (20,30)

2.射線判別法

 

根據對多邊形的瞭解，我們可以得出如下結論：

如果一個點在多邊形內部，任意角度做射線肯定會與多邊形要麼有一個交點，要麼有與多邊形邊界線重疊。

如果一個點在多邊形外部，任意角度做射線要麼與多邊形有一個交點，要麼有兩個交點，要麼沒有交點，要麼有與多邊形邊界線重疊。

利用上面的結論，我們只要判斷這個點與多邊形的交點個數，就可以判斷出點與多邊形的位置關係了。

首先羅列下注意事項:

l 射線跟多邊形的邊界線重疊的情況

l 區別內部點和外部點的射線在有一個交點時的情況

對於第一個注意事項，可以將射線角度設為零度，這樣子只需要判斷兩個相鄰頂點的Y值是否相等即可。然後再判斷這個點的X值方位。

對於第二個注意事項，網上許多文章都說到做射線以後交點為奇數則表示在多邊形內部，這是一個錯誤的觀點，不僅對於凹多邊形不成立，對於凸多邊形也不成立。

例如：從外部點做射線剛好經過一頂點，這樣子交點個數就為奇數，但是該點卻不在多邊形內部。

至於要如何區分這兩種情況呢，我能想到了一個不完美的方法，外部點的射線跟多邊形有一個交點的時候，該交點肯定為頂點，如果該射線上移一位或者下移一位，要麼變成有兩個交點要麼沒有交點。當然為了安全起見，這裡把射線盡量往相鄰點中心移動。這樣子就能夠判斷出是外部點的射線跟多邊形有一個交點。

不過這個方法並不完美，因為有了移位操作，可能會把內部點移動出外部。而且如果判斷點在(60,30)位置，判斷的時候先遇到(20,30)，然後移位操作，就判斷就出錯了。

為瞭解決這些問題，我在起初先掃描一次判斷點是否在頂點上雖然影響了一點效率，而且當判定點距離多邊形一個單位時，判斷可能會有誤。不過只要不是需要高精度的話，這個方法還是很有效。

代碼如下：

 

bool InsidePolygon1( POINTD *polygon,int N,POINTD pt ){    int i,j;    bool inside,redo;    redo = true;    for (i = 0;i < N;++i)    {        if (polygon[i].x == pt.x &&    // 是否在頂點上            polygon[i].y == pt.y )        {            redo = false;            inside = true;            break;        }    }    while (redo)    {        redo = false;        inside = false;        for (i = 0,j = N - 1;i < N;j = i++)         {            if ( (polygon[i].y < pt.y && pt.y < polygon[j].y) ||                 (polygon[j].y < pt.y && pt.y < polygon[i].y) )             {                if (pt.x <= polygon[i].x || pt.x <= polygon[j].x)                 {                    double _x = (pt.y-polygon[i].y)*(polygon[j].x-polygon[i].x)/(polygon[j].y-polygon[i].y)+polygon[i].x;                    if (pt.x < _x)          // 線上的左側                        inside = !inside;                    else if (pt.x == _x)    // 線上上                    {                        inside = true;                        break;                    }                }            }            else if ( pt.y == polygon[i].y)             {                if (pt.x < polygon[i].x)    // 交點在頂點上                {                    polygon[i].y > polygon[j].y ? --pt.y : ++pt.y;                    redo = true;                    break;                }            }            else if ( polygon[i].y ==  polygon[j].y && // 在水平的邊界線上                pt.y == polygon[i].y &&                ( (polygon[i].x < pt.x && pt.x < polygon[j].x) ||                 (polygon[j].x < pt.x && pt.x < polygon[i].x) ) )            {                inside = true;                break;            }        }    }    return inside;}

 

3.角度和判別法

 

  角度和判別法就是判定點與多邊形所有相鄰頂點組成的角的角度和來判斷點與多邊形的位置關係。

如果點在多邊形內部，只要該點不在邊界線或者頂點上，則角度和為三百六十度。

如果在邊界線上，則角度和為一百八十度。

如果點在多邊形外部，則角度和無法達到三百六十度。但是角度和可能會達到一百八十度，比如判斷點在(60,10)。

代碼如下：

 

// 根據需要不判斷頂點bool IsPointInLine( const POINTD &pt,const POINTD &pt1,const POINTD &pt2 ){    bool inside = false;    if (pt.y == pt1.y &&        pt1.y == pt2.y &&        ((pt1.x < pt.x && pt.x < pt2.x) ||         (pt2.x < pt.x && pt.x < pt1.x)) )    {        inside = true;    }    else if (pt.x == pt1.x &&        pt1.x == pt2.x &&        ((pt1.y < pt.y && pt.y < pt2.y) ||         (pt2.y < pt.y && pt.y < pt1.y)) )    {        inside = true;    }    else if ( ((pt1.y < pt.y && pt.y < pt2.y) ||         (pt2.y < pt.y && pt.y < pt1.y)) &&        ((pt1.x < pt.x && pt.x < pt2.x) ||         (pt2.x < pt.x && pt.x < pt1.x)) )    {        if (0 == (pt.y-pt1.y)/(pt2.y-pt1.y)-(pt.x - pt1.x) / (pt2.x-pt1.x))        {            inside = true;        }    }    return inside;}bool InsidePolygon2( POINTD *polygon,int N,POINTD p ){    int i,j;    double angle = 0;    bool inside = false;    for (i = 0,j = N - 1;i < N;j = i++)     {        if (polygon[i].x == p.x &&    // 是否在頂點上            polygon[i].y == p.y)        {            inside = true;            break;        }        else if (IsPointInLine(p,polygon[i],polygon[j]))    // 是否在邊界線上        {            inside = true;            break;        }        double x1,y1,x2,y2;        x1 = polygon[i].x - p.x;        y1 = polygon[i].y - p.y;        x2 = polygon[j].x - p.x;        y2 = polygon[j].y - p.y;        double radian = atan2(y1,x1) - atan2(y2,x2);        radian = abs(radian);        if (radian > M_PI) radian = 2* M_PI - radian;        angle += radian;            // 計算角度和    }    if ( fabs(6.28318530717958647692 - angle) < 1e-7 )        inside = true;    return inside;}

 

有的人管角度和判別法叫做轉角法，還有一個類似的方法叫弧長法。

http://haibarali.blog.163.com/blog/static/23474013200722813356858/ 

4.面積和判別法

 

根據角度和判別法，我們可以繼續演化，可以得出如下結論：

如果一個點在多邊形內部，該點與多邊形所有相鄰頂點組成的三角形面積和為多邊形面積。反之不成立。

求三角形面積：

已知三角形A（x1,y1）B（x2,y2）C（x3,y3），則面積公式為：

                 |x1  x2  x3|   

S(A,B,C) = |y1  y2  y3| * 0.5 (當三點為逆時針時為正，順時針則為負的)   

                 |1   1   1 |  

= ( x1*y2 - x1*y3 - x2*y1 + x2*y3 + x3*y1 - x3*y2 ) * 0.5

求多邊形面積：

根據上面求三角形的方法，

已知多邊形A1A2A3...An（順或逆時針都可以），設平面上有任意的一點P，則面積公式為：   

S(A1,A2,A3...An)   

     =  S(P，A1，A2)+ S(P，A2，A3)+...+S(P,An,A1)

    

既然P是可以任取，為了方便用(0,0)好了。

代碼如下：

bool InsidePolygon3( POINTD *polygon,int N,POINTD pt ){    int i,j;    bool inside = false;    double polygon_area = 0;    double trigon_area = 0;    for (i = 0,j = N - 1;i < N;j = i++)     {        polygon_area += polygon[i].x * polygon[j].y - polygon[j].x * polygon[i].y;        trigon_area += abs(            pt.x * polygon[i].y -            pt.x * polygon[j].y -            polygon[i].x * pt.y +            polygon[i].x * polygon[j].y +            polygon[j].x * pt.y -            polygon[j].x * polygon[i].y            );    }    trigon_area *= 0.5;    polygon_area = abs(polygon_area * 0.5);    if ( fabs(trigon_area - polygon_area) < 1e-7 )        inside = true;    return inside;}

多邊形面積計算

http://www.cppblog.com/zyzx/archive/2009/04/27/81241.html 

行列式如何展開

http://www.tongji.edu.cn/~math/xxds/kcja/kcja_b/1-6.htm 

5.點線判別法

經網友glshader提示，添加了這個方法。

如果判斷點在所有邊界線的同側，就能判定該點在多邊形內部。

判斷方法就是判斷兩條同起點射線斜率差。

代碼如下：

 

bool InsidePolygon4( POINTD *polygon,int N,POINTD p ){    int i,j;    bool inside = false;    int count1 = 0;    int count2 = 0;    for (i = 0,j = N - 1;i < N;j = i++)     {        double value = (p.x - polygon[j].x) * (polygon[i].y - polygon[j].y) - (p.y - polygon[j].y) * (polygon[i].x - polygon[j].x);        if (value > 0)            ++count1;        else if (value < 0)            ++count2;    }    if (0 == count1 ||        0 == count2)    {        inside = true;    }    return inside;}

 

6.總結

    總結了四種判別法，雖然很多原理都很簡單，但是實際操作起來，發現很多細節還是值得推敲的。

我在再來個原創java版本的，不過沒有射線的那個方法準確！

/**   * 判定一個經緯度點與地區的關係   * @param point 需要判斷的經緯度點   * @return boolean 如果點在地區內，true；如果在地區外，false   */  public boolean aboutZone(Point point)  {    if(vPoint==null||vPoint.size()==0) return false;    int I = 0;    //求出最大值, 得到一個落在矩形外的點, 由該點與point組成的線段, 與多邊形求交, 交點為奇數則為    // 在多邊形內    I = 1;    int odd = 0;    long ytemp;    Point pt0,pt1, pt2, pt3, pt4;    while( I < vPoint.size()-2 )    {      pt0 = (Point)vPoint.get(I - 1);      pt1 = (Point)vPoint.get(I);      pt2 = (Point)vPoint.get(I+1);      pt3 = (Point)vPoint.get(I+2);      if ( ( (point.lon > pt1.lon) &&             (point.lon < pt2.lon)) ||             ( (point.lon < pt1.lon) &&             (point.lon > pt2.lon))             ) {        ytemp = pt1.lat +                (point.lon - pt1.lon) *  (pt2.lat - pt1.lat) /                (pt2.lon - pt1.lon);        if (  point.lat > ytemp)          odd++;        else if (  point.lat == ytemp)          return true;      }      else if (pt1.lon == point.lon) { // begin of special cases (on sides)        if (point.lat == pt1.lat)          return true;        if (pt1.lon == pt2.lon) {          if ( (point.lat - pt1.lat) * (point.lat - pt2.lat) <= 0)            return true;          if ( (point.lon - pt0.lon) * (point.lon - pt3.lon) < 0) {            if (point.lat > pt1.lat)              odd++;          }        }        else          if ( ( (point.lon - pt0.lon) * (point.lon - pt2.lon) <=0 ) &&               (point.lat > pt1.lat)) {            odd++;          }      }      I ++;    }    pt0 = (Point)vPoint.get(0);    pt1 = (Point)vPoint.get(1);    pt2 = (Point)vPoint.get(2);    pt3 = (Point)vPoint.get(vPoint.size() - 1);    pt4 = (Point)vPoint.get(vPoint.size() - 2);    if ( ( (point.lon > pt0.lon) &&           (point.lon < pt1.lon)) ||           ( (point.lon < pt0.lon) &&           (point.lon > pt1.lon))) {      ytemp = pt0.lat +              (point.lon - pt0.lon) * (pt1.lat - pt0.lat) /(pt1.lon - pt0.lon);      if (point.lat > ytemp)        odd++;      else if (point.lat == ytemp)        return true;    }    else if (pt0.lon == point.lon) { // begin of special cases (on sides)      if (point.lat == pt0.lat)        return true;      if (pt0.lon == pt1.lon) {        if ( (point.lat - pt0.lat) * (point.lat - pt1.lat) <= 0.1f)          return true;        if ( (point.lon - pt3.lon) * (point.lon - pt2.lon) < 0.f) {          if (point.lat > pt0.lat) {            odd++;          }        }      }      else if ( ( (point.lon - pt3.lon) * (point.lon - pt1.lon) < 0.1f) &&                (point.lat > pt0.lat)) {        odd++;      }    }    if ( ( (point.lon > pt3.lon) && (point.lon < pt0.lon)) || ( (point.lon < pt3.lon) &&           (point.lon > pt0.lon))) {      ytemp = pt3.lat +              (point.lon - pt3.lon) * (pt0.lat - pt3.lat) /(pt0.lon - pt3.lon);      if (point.lat > ytemp)        odd++;      else if (  point.lat == ytemp)        return true;    }    else if (pt3.lon == point.lon) {      if (point.lat == pt3.lat)        return true;      if (pt3.lon == pt0.lon) {        if ( (point.lat - pt3.lat) * (point.lat - pt0.lat) <= 0)          return true;        if ( (point.lon - pt4.lon) * (point.lon - pt1.lon) <= 0 ){          if (point.lat > pt3.lat) {            odd++;          }        }      }      else if ( ( (point.lon - pt4.lon) * (point.lon - pt0.lon) <= 0) &&                (point.lat > pt3.lat)) {        odd++;      }    }    if ( ( (odd >> 1) << 1) == odd)      return false;    return true;  } 

package com.encryptionalgorithm;public class Point {private long longitude;private long latitude ;private double x;private double y;public void setX(double x){this.x = x;}public void setY(double y){this.y = y;}public double getX(){return x;}public double getY(){return y;}public void setLongitude(long longitude){this.longitude = longitude; }public void setLatitude(long latitude){this.latitude = latitude;}public long getLongitude(){return longitude;}public long getLatitude(){return latitude;}}