点与多边形关系（改进射线法）

在GIS软件开发中，经常要用到一些几何的算法，比如三角网构建，多边形的剖分，点，线，面之间的关系。而点与多边形关系的判断是一项非常重要的基础工作。

在点与多边形关系的判断中，经常用到的方法是射线法和夹角和方法，其中射线法能够针对带岛的多边形进行判断，而夹角和方法就显得无能为力。

射线法的基本思想是：从待判断的点向某一个方向引射线，计算和多边形交点的个数，如果个数是偶数或者0则点在多边形外，如果是奇数，则在多边形内。这个只是最基本的判别情况，还有一些复杂的情况需要特殊处理：

（射线经过顶点）：当射线经过顶点时，判断就会出现异常情况，现在规定，线段的两个端点，相对于另一个端点在上面的顶点称为上端点，下面是下端点，如果经过下端点，则认为边和射线不相交。

（点在边上）：这种情况也不能用交点个数的奇偶性来判断了，要快速地判断这个点是否在边上。

射线法改进：传统的射线法一开始就直接计算点和多边形的交点个数，这样的话，会花费大量的时间来作拓扑关系的判断。改进的算法是首先利用多边形的最小外接矩形迅速排出掉不在MBR内的点，然后利用交点个数的奇偶性判断：

下面的函数是射线和边关系以及交点个数判断：

//射线和线段的关系 :相交返回1，不相交返回0，射线起点在线段上返回-1
int IsIntersectAnt(double x,double y,double X1,double Y1,double X2,double Y2)
{
	//计算线段的最小和最大坐标值
	double minX,maxX,minY,maxY;
	minX = X1;
	maxX = X2;
	if (minX > maxX)
	{
		minX = X2;
		maxX = X1;
	}
	minY = Y1;
	maxY = Y2;
	if (minY > maxY)
	{
		minY = Y2;
		maxY = Y1;
	}

	//射线与边无交点的快速判断
	if (y<minY || y>maxY || x<minX)
	{
		return 0;
	}

	//如果是水平线段，在线段上返回-1，否则返回0
	if (fabs(maxY - minY) < eps)
	{
		return (x >= minX && x <= maxX)? (-1):0;
	}

	//计算射线与边所在直线的交点的横坐标
	double x0 = X1 + (double)(y - Y1)*(X2 - X1)/(Y2 - Y1);
	
	//交点在射线右侧，则不相交
	if (x0 > x)
	{
		return 0;
	}
	//交点和射线起点相同
	if (fabs(x-x0)< eps)
	{
		return -1;
	}
	//穿过下端点也不计数
	if (fabs(y-minY) < eps)
	{
		return 0;
	}
	return 1;

}

上面的是计算射线和一条边的情况，对于一个多边形所以只要逐个判断射线和它的边就可以了

int MyPolygon::PointInPolygon(const MyPoint& poPoint)
{
	//如果点不在多边形的最小外接矩形中，则一定不在多边形内
	MyEnvelope env;
	GetEnvelope(&env);
	if (!env.IsPointInRect(poPoint.GetX(),poPoint.GetY()))
	{
		return 0;
	}

	//计算该点向左方向的射线与各个边的交点个数
	int nCount = 0;
	double X = poPoint.GetX();
	double Y = poPoint.GetY();
	int nFlag = 0;
	MyLineString *exteriorRing = GetExteriorRing();
	exteriorRing->CloseRings();
	for (int i = 0; i < exteriorRing->GetNumPoint()-1; i ++)
	{
		nFlag = IsIntersectAnt(X,Y,exteriorRing->GetX(i),exteriorRing->GetY(i),
			exteriorRing->GetX(i+1),exteriorRing->GetY(i+1));
		if (nCount < 0) 
		{
			return 2;	//点在边上
		}
		nCount += nFlag;
	}

	if (nCount % 2 == 1)   //点在多边形内
	{
		return 1;
	}
	else
		return 0;
}

程序经过测试是没有问题的，欢迎大家批评指正！