#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;
}