ptinpolygon.cpp

测试点是否落在多边形内

#include "shapefil.h"
#include "PtInPolygon.h"

int  PointInPolygon(const SHPObject *polygonObject, double x, double y)
{
//////////////////////////////////////////////
//											//
//		点与多边形位置关系的分析算法		//
//											//
//////////////////////////////////////////////

//基本思想（射线算法）
//计算从测试点发出的向上的射线与多边形边的交点总数，
//如果交点总数为奇数，说明点在多边形内，返回值为2，
//否则说明点在多边形外，返回值为0；如果点在多边形边
//上，返回值为1。

//特殊情况处理
//1、如果射线穿过多边形的顶点，需判断经过该顶点的两
//条边与射线的位置关系，如果两条边在射线的同一侧，此
//时算两个交点；如果在不同侧，此时只算一个交点
//2、如果多边形的边是垂直线段、位于测试点上方，且与
//测试点的坐标重合，此时把该边视作顶点来处理。

	long nVertices = polygonObject->nVertices;
	if(nVertices < 3)
		return 0;		//Not a valid polygon, return 0
	UINT returnvalue = -1;
	int intersectPtNum = 0;
	BOOL bFirst = TRUE;
	int startPtIndex = 0;
	Point pt = Point(x, y);

	for(int i = 0, p=0; i < nVertices-1 && p < polygonObject->nParts; i++)
	{
		if(bFirst)
		{
			startPtIndex = i;
			bFirst = FALSE;
		}
		if((p < polygonObject->nParts-1 && i == polygonObject->panPartStart[p+1]))
		{
			bFirst = TRUE;
			p++;
			continue;
		}
//		if(startPtIndex != i && pt1 == Point(polygonObject->padfX[startPtIndex], polygonObject->padfY[startPtIndex]))
//		//如果多边形外边界或内岛边界已经分析完成，转入下一岛的分析
//		{
//			bFirst = TRUE;
//			continue;
//		}
		Point pt1, pt2;
		pt1 = Point(polygonObject->padfX[i], polygonObject->padfY[i]);
		pt2 = Point(polygonObject->padfX[i+1], polygonObject->padfY[i+1]);
		//获取多边形边的最大和最小x取值范围，如果测试点的x坐标位于该范围之外
		//则跳到另一边继续测试
		double minx, maxx;	
		minx = pt1.x > pt2.x ? pt2.x : pt1.x;
		maxx = pt1.x < pt2.x ? pt2.x : pt1.x;
		if(minx - pt.x > Global_Tolerance || pt.x - maxx > Global_Tolerance)
			continue;
		if(fabs(pt1.x - pt2.x) < Global_Tolerance)
		//如果多边形的边是垂直的，并且其x值与测试点相等
		//则判断点是否在这条边上，如果在边上，则返回值等于1，完成测试
		{
			if(fabs(pt.x - pt1.x) < Global_Tolerance)
			//如果测试点的x坐标等于垂直边的x坐标（多余之举，完整算法应包含）
			{
				//获取垂直边的最大和最小y取值范围
				//如果测试点的y值在这个范围之内，
				//说明点在直线上，返回值等于1
				//否则说明点不在直线上
				double miny, maxy;
				miny = pt1.y > pt2.y ? pt2.y : pt1.y;
				maxy = pt1.y < pt2.y ? pt2.y : pt1.y;
				if(pt.y - maxy < Global_Tolerance && miny - pt.y < Global_Tolerance) 
				{
					returnvalue = 1;
					break;
				}
				else
					continue;
			}
			else
				continue;
		}
		//如果多边形的边不垂直，计算其斜率
		//并计算穿过测试点的垂直直线与该边的交点y坐标
		double k, b, y;
		k = (pt1.y - pt2.y)/(pt1.x - pt2.x);
		b = pt1.y - pt1.x*k;
		y = pt.x*k + b;
		if(fabs(y - pt.y) < Global_Tolerance)
		//如果y坐标与测试点相同，说明测试点在这条边上；
		{
			returnvalue = 1;
			break;
		}
		if(y > pt.y)
		//如果y坐标大于测试点的y坐标，说明从测试点发出的向上的射线与该边存在交点
		{
			if(fabs(pt.x - pt1.x) < Global_Tolerance)
			//如果射线穿过多边形边的起点
			{
				Point pt11, pt22;
				int j = i-1;
				while(TRUE)
				//找出该顶点之前第一个非垂直的边的顶点
				{
					if(j < 0)
						j = nVertices-2;
					if(fabs(polygonObject->padfX[j] - pt.x) > Global_Tolerance)
					{
						pt11 = Point(polygonObject->padfX[j], polygonObject->padfY[j]); break;
					}
					j--;
				}
				pt22 = Point(polygonObject->padfX[i+1], polygonObject->padfY[i+1]);

				if((pt.x - pt11.x > Global_Tolerance &&
					pt.x - pt22.x > Global_Tolerance) ||
					(pt11.x - pt.x > Global_Tolerance &&
					pt22.x - pt.x > Global_Tolerance))
					intersectPtNum += 2;
				else
					intersectPtNum++;
			}
			else if(fabs(pt.x - pt2.x) > Global_Tolerance)
				intersectPtNum++;
		}
	}
	if(returnvalue == -1)
		returnvalue = (intersectPtNum%2 == 0) ? 0 : 2;
	return returnvalue;
}