判断一个点是否在多边形内

#转载自: http://blog.csdn.net/u011722133/article/details/52813374

在GIS（地理信息管理系统)/PCL(点云库)中，判断一个坐标是否在多边形内部是个经常要遇到的问题。乍听起来还挺复杂。根据W. Randolph Franklin 提出的PNPoly算法，只需区区几行代码就解决了这个问题

假设多边形的坐标存放在一个数组里，首先我们需要取得该数组在横坐标和纵坐标的最大值和最小值，根据这四个点算出一个四边型，首先判断目标坐标点是否在这个四边型之内，如果在这个四边型之外，那可以跳过后面较为复杂的计算;

if (p.x < minX || p.x > maxX || p.y < minY || p.y > maxY) 
{
　　//如果难以通过筛选，直接返回false
}

接下来是核心算法部分

int pnpoly (int nvert, float *vertx, float *verty, float testx, float testy) 
{
    int i, j;
    bool c = 0;
    for (i = 0, j = nvert-1; i < nvert; j = i++) 
    {
　　　　 bool cond_1 =(verty[i]>testy) != (verty[j]>testy);
　　　　 bool cond_2 =(testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]));

　　　　if ( cond_1 && cond_2}
　　　　{
　　　　　　c = !c; 
        }
    }
    return c;
}

首先，参数nvert 代表多边形有几个点。浮点数testx, testy代表待测试点的横坐标和纵坐标，*vertx,*verty分别指向储存多边形横纵坐标数组的首地址。
我们注意到，每次计算都涉及到相邻的两个点和待测试点，然后考虑两个问题：
1. 被测试点的纵坐标testy是否在本次循环所测试的两个相邻点纵坐标范围之内？即verty[i] <testy < verty[j] 或者 verty[j] <testy < verty[i]
2. 待测点test是否在i,j两点之间的连线之下？利用i,j两点间的斜率公式，然后每次这两个条件同时满足的时候我们把返回的布尔量取反。

这个表达式的意思是说，随便画个多边形，随便定一个点，然后通过这个点水平划一条射线，先数数看这条射线和多边形的边相交几次，（或者说先排除那些不相交的边，第一个判断条件），然后再数这条射
线穿越多边形的次数是否为奇数，如果是奇数，那么该点在多边形内，如果是偶数，则在多边形外。

转载阿凡卢，转载地址：http://www.cnblogs.com/luxiaoxun/p/3722358.html

 理论参照 

判断点是否在多边形内部
如何判断一个点是否在多边形内部？
(1) 面积和判别法：判断目标点与多边形的每条边组成的三角形面积和是否等于该多边形，相等则在多边形内部。
(2) 夹角和判别法：判断目标点与所有边的夹角和是否为360度，为360度则在多边形内部。
(3) 引射线法：从目标点出发引一条射线，看这条射线和多边形所有边的交点数目。如果有奇数个交点，则说明在内部，如果有偶数个交点，则说明在外部。
具体做法：将测试点的Y坐标与多边形的每一个点进行比较，会得到一个测试点所在的行与多边形边的交点的列表。在下图的这个例子中有8条边与测试点所在的行相交，而有6条边没有相交。如果测试点
的两边点的个数都是奇数个则该测试点在多边形内，否则在多边形外。在这个例子中测试点的左边有5个交点，右边有三个交点，它们都是奇数，所以点在多边形内。

算法图解

关于这个算法的具体的更多图形例子：http://alienryderflex.com/polygon/

参考代码

int pnpoly(int nvert, float *vertx, float *verty, float testx, float testy)
{
  int i, j, c = 0;
  for (i = 0, j = nvert-1; i < nvert; j = i++) 
  {
    if ( ((verty[i]>testy) != (verty[j]>testy)) &&
     (testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) )
       c = !c;
  }
  return c;
}

一个来自polygon内部的实现

public bool IsInside(PointLatLng p)
      {
         int count = Points.Count;

         if(count < 3)
         {
            return false;
         }

         bool result = false;

         for(int i = 0, j = count - 1; i < count; i++)
         {
            var p1 = Points[i];
            var p2 = Points[j];

            if(p1.Lat < p.Lat && p2.Lat >= p.Lat || p2.Lat < p.Lat && p1.Lat >= p.Lat)
            {
               if(p1.Lng + (p.Lat - p1.Lat) / (p2.Lat - p1.Lat) * (p2.Lng - p1.Lng) < p.Lng)
               {
                  result = !result;
               }
            }
            j = i;
         }
         return result;
      }

特殊情况：要检测的点在多变形的一条边上，射线法判断的结果是不确定的，需要特殊处理（If the test point is on the border of the polygon, this algorithm will deliver unpredictable results)
计算一个多边形的面积（area of a polygon）：

private static double SignedPolygonArea(List<PointLatLng> points)
        {
            // Add the first point to the end.
            int pointsCount = points.Count;
            PointLatLng[] pts = new PointLatLng[pointsCount + 1];
            points.CopyTo(pts, 0);
            pts[pointsCount] = points[0];

            for (int i = 0; i < pointsCount + 1; ++i)
            {
                pts[i].Lat = pts[i].Lat * (System.Math.PI * 6378137 / 180);
                pts[i].Lng = pts[i].Lng * (System.Math.PI * 6378137 / 180);
            }

            // Get the areas.
            double area = 0;
            for (int i = 0; i < pointsCount; i++)
            {
                area += (pts[i + 1].Lat - pts[i].Lat) * (pts[i + 1].Lng + pts[i].Lng) / 2;
            }

            // Return the result.
            return area;
        }

        /// <summary>
        /// Get the area of a polygon
        /// </summary>
        /// <param name="points"></param>
        /// <returns></returns>
        public static double GetPolygonArea(List<PointLatLng> points)
        {
            // Return the absolute value of the signed area.
            // The signed area is negative if the polygon is oriented clockwise.
            return Math.Abs(SignedPolygonArea(points));
        }