GCS-GISControlDlg-for-981A-and-981C/TopologicalAnalysis.cpp

#include "stdafx.h"
#include "topologicalanalysis.h"


TopologicalAnalysis::TopologicalAnalysis(void)
{
}


TopologicalAnalysis::~TopologicalAnalysis(void)
{
}


bool TopologicalAnalysis::isPointInLine(double* point, double* startPoint, double* endPoint,float tolerance)
{
	double AX = startPoint[0];
	double AY = startPoint[1];
	double BX = endPoint[0];
	double BY = endPoint[1];
	double PX = point[0];
	double PY = point[1];

	double dx_AB = AX - BX;
	double dy_AB = AY - BY;
	double dx_PA = PX - AX;
	double dy_PA = PY - AY;
	double dx_PB = PX - BX;
	double dy_PB = PY - BY;

	double AB = sqrt(dx_AB*dx_AB + dy_AB*dy_AB);
	double PA = sqrt(dx_PA*dx_PA + dy_PA*dy_PA);
	double PB = sqrt(dx_PB*dx_PB + dy_PB*dy_PB);
	double rate = abs(PA + PB - AB) / AB;
	if (rate < tolerance)
	{
		return true;
	}
	else
	{
		return false;
	}
}

//<2F>ж<EFBFBD><D0B6><EFBFBD><EFBFBD>Ƿ<EFBFBD><C7B7><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
int TopologicalAnalysis::isPointInPolyLine(double* point, vector<double>& lineX,vector<double>& lineY,float tolerance)
{
	int lineNum = lineX.size();
	double startPoint[2],endPoint[2];
	for (int i=0;i<lineNum-1;i++)
	{
		startPoint[0] = lineX.at(i);
		startPoint[1] = lineY.at(i);
		endPoint[0] = lineX.at(i+1);
		endPoint[1] = lineY.at(i+1);
		bool b_in = isPointInLine(point,startPoint,endPoint,tolerance);
		if(b_in)
		{
			return (i+1);
		}
	}
	//<2F>ж<EFBFBD><D0B6><EFBFBD>β<EFBFBD><CEB2><EFBFBD>߶<EFBFBD>
	startPoint[0] = lineX.at(lineNum-1);
	startPoint[1] = lineY.at(lineNum-1);
	endPoint[0] = lineX.at(0);
	endPoint[1] = lineY.at(0);
	bool b_end = isPointInLine(point,startPoint,endPoint,tolerance);
	if (b_end)
	{
		return lineNum;
	}
	else
	{
		return 0;
	}
}

//<2F><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>߹<EFBFBD><DFB9><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֱ<EFBFBD>ߵĽ<DFB5><C4BD><EFBFBD>
bool TopologicalAnalysis::GetPointToLineVerticalCross(double* linePt1,double* linePt2,double* pt,double* crossPt)
{
	//<2F><>ֱ<EFBFBD><D6B1>
	if (linePt1[0] == linePt2[0])
	{
		crossPt[0] = linePt1[0];
		crossPt[1] = pt[1];
		return true;
	}
	//ˮƽ<CBAE><C6BD>
	if (linePt1[1] == linePt2[1])
	{
		crossPt[0] = pt[0];
		crossPt[1] = linePt1[1];
		return true;
	}
	float A = (linePt1[1]- linePt2[1]) * 1.0 / (linePt1[0]- linePt2[0]);
	float B = (linePt1[1] - A * linePt1[0]);
	float m = pt[0] + A * pt[1];

	/// <20><><EFBFBD><EFBFBD>ֱ<EFBFBD>߽<EFBFBD><DFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
	crossPt[0] = (m - A * B) * 1.0f / (A * A + 1);
	crossPt[1] = A * crossPt[0] + B;
	return true;
}

//<2F>жϵ<D0B6><CFB5>Ƿ<EFBFBD><C7B7><EFBFBD>(<28><><EFBFBD><EFBFBD>)<29><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
//polygon: <20><>β<EFBFBD><CEB2>ͬ<EFBFBD><CDAC>Cpoint1<74>б<EFBFBD>
bool TopologicalAnalysis::isPointInPolygon(CPoint1 point, vector<CPoint1> polygon){
	if (polygon.size()<=3) return false;	// һ<><D2BB><EFBFBD><EFBFBD>Ч<EFBFBD><D0A7><EFBFBD><EFBFBD><EFBFBD>ζ<EFBFBD><CEB6><EFBFBD><EFBFBD><EFBFBD>Ӧ<EFBFBD><D3A6><EFBFBD><EFBFBD>3

	int LineNum = polygon.size();
	CPoint1 leftP = point;
	CPoint1 rightP;
	rightP.SetX(getMaxX(polygon) + 1);
	rightP.SetY(point.GetY());
	int count = 0, yPrev = polygon[LineNum - 2].GetY();
	CPoint1 v1, v2;
	v1 = polygon[LineNum - 1];
	for (int i = 0; i < LineNum; i++)
	{
		v2 = polygon[i];

		if (isPointInLine(leftP, v1, v2))
			return true;

		if (v1.GetY() != v2.GetY())
		{
			if (isLineIntersect(v1, v2, leftP, rightP))
			{
				if (isPointInLine(v1, leftP, rightP))
				{
					if (v1.GetY()<v2.GetY()) { if (v1.GetY()>yPrev)  count++; }
					else { if (v1.GetY() < yPrev) count++; }
				}
				else if (!isPointInLine(v2, leftP, rightP))
				{
					count++;
				}
			}
		}
		yPrev = v1.GetY();

		v1 = v2;
	}
	return   (count % 2 == 1);
}

double TopologicalAnalysis::getMaxX(vector<CPoint1> points){
	if (points.size()==0)
		return -1;
	else if(points.size()==1)
		return points[0].GetX();
	else{
		double maxx = points[0].GetX();
		for (unsigned i=1; i<points.size();i++)
		{
			if(points[i].GetX()>maxx){
				maxx = points[i].GetX();
			}
		}
		return maxx;
	}
}

bool  TopologicalAnalysis::isPointInLine(CPoint1 point, CPoint1 startPoint, CPoint1 endPoint)
{
	long AX = startPoint.GetX();
	long AY = startPoint.GetY();
	long BX = endPoint.GetX();
	long BY = endPoint.GetY();
	long PX = point.GetX();
	long PY = point.GetY();

	double dx_AB = AX - BX;
	double dy_AB = AY - BY;
	double dx_PA = PX - AX;
	double dy_PA = PY - AY;
	double dx_PB = PX - BX;
	double dy_PB = PY - BY;

	double AB = sqrt(dx_AB*dx_AB + dy_AB*dy_AB);
	double PA = sqrt(dx_PA*dx_PA + dy_PA*dy_PA);
	double PB = sqrt(dx_PB*dx_PB + dy_PB*dy_PB);
	double rate = abs(PA + PB - AB) / AB;
	if (rate < 0.001)
	{
		return true;
	}
	else
	{
		return false;
	}
}
// <20><><EFBFBD><EFBFBD>
double mult(CPoint1 a, CPoint1 b, CPoint1 c)  
{  
	return (a.GetX()-c.GetX())*(b.GetY()-c.GetY())-(b.GetX()-c.GetX())*(a.GetY()-c.GetY());  
}

//<2F>ж<EFBFBD><D0B6><EFBFBD><EFBFBD><EFBFBD>ֱ<EFBFBD>߶<EFBFBD><DFB6>Ƿ<EFBFBD><C7B7>ཻ
bool TopologicalAnalysis::isLineIntersect(CPoint1 line1Start, CPoint1 line1End, CPoint1 line2Start, CPoint1 line2End){
	double l1sx = line1Start.GetX();
	double l1sy = line1Start.GetY();
	double l1ex = line1End.GetX();
	double l1ey = line1End.GetY();
	double l2sx = line2Start.GetX();
	double l2sy = line2Start.GetY();
	double l2ex = line2End.GetX(); 
	double l2ey = line2End.GetY();
	if ( max(l1sx, l1ex)<min(l2sx, l2ex) )  
	{  
		return false;  
	}  
	if ( max(l1sy, l1ey)<min(l2sy,l2ey) )  
	{  
		return false;  
	}  
	if ( max(l2sx, l2ex)<min(l1sx, l1ex) )  
	{  
		return false;  
	}  
	if ( max(l2sy,l2ey)<min(l1sy, l1ey) )  
	{  
		return false;  
	}  
	if ( mult(line2Start, line1End, line1Start)*mult(line1End, line2End, line1Start)<0 )  
	{  
		return false;  
	}  
	if ( mult(line1Start, line2End, line2Start)*mult(line2End, line1End, line2Start)<0 )  
	{  
		return false;  
	}  
	return true;  
}


bool TopologicalAnalysis::GetBoundingBoxVertices(const vector<double>& polygonX,const vector<double>& polygonY,vector<double>& rectangleX,vector<double>& rectangleY)
{
	if (polygonX.size()<3 || polygonY.size()<3) return false;

	//<2F><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Сx<D0A1><78>y
	double minX = polygonX[0];
	double minY = polygonY[0];
	double maxX = polygonX[0];
	double maxY = polygonY[0];

	for (int i=0;i<polygonY.size();++i) {
		if (polygonX[i] < minX) minX = polygonX[i];
		if (polygonY[i] < minY) minY = polygonY[i];
		if (polygonX[i] > maxX) maxX = polygonX[i];
		if (polygonY[i] > maxY) maxY = polygonY[i];
	}

	rectangleX.push_back(minX);
	rectangleX.push_back(maxX);
	rectangleX.push_back(maxX);
	rectangleX.push_back(minX);

	rectangleY.push_back(maxY);
	rectangleY.push_back(maxY);
	rectangleY.push_back(minY);
	rectangleY.push_back(minY);
}

//<2F><><EFBFBD><EFBFBD>ֱ<EFBFBD><D6B1><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>εĽ<CEB5><C4BD><EFBFBD>
int TopologicalAnalysis::linePolygonIntersections(const Point2D& linePt1,const Point2D& linePt2,const vector<double>& polygonX,const vector<double>& polygonY, vector<Point2D>& result)
{

	Line2D line1;
	line1.p1 = linePt1;
	line1.p2 = linePt2;
	line1.is_seg = false;
	Point2D intersectionPt;
	vector<Point2D> resPoints;
	for (int i=0;i<polygonX.size()-1;++i)
	{
		Point2D pt1,pt2;
		pt1.x = polygonX[i];
		pt1.y = polygonY[i];
		pt2.x = polygonX[i+1];
		pt2.y = polygonY[i+1];
		Line2D line2;
		line2.p1 = pt1;
		line2.p2 = pt2;
		line2.is_seg = false;

/*
		int type = isLineIntersecting(line1,line2,intersectionPt);
		if (type == 1 || type == 2) //<2F>ཻ<EFBFBD><E0BDBB><EFBFBD>غ<EFBFBD>
		{
			resPoints.push_back(intersectionPt);
		}*/

/*
		if (isSegIntersect(line1,line2,intersectionPt))
		{
			line2.is_seg = true;
			if (isponl(intersectionPt,line2))
			{
				resPoints.push_back(intersectionPt);
			}	
			//resPoints.push_back(intersectionPt);
		}*/

		
		if (findIntersection(linePt1.x,linePt1.y,linePt2.x,linePt2.y,pt1.x,pt1.y,pt2.x,pt2.y,intersectionPt))
		{
			resPoints.push_back(intersectionPt);
		}
	}

	if (resPoints.size()>1)
	{
		//<2F><>ȡ<EFBFBD><C8A1><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ҵ<EFBFBD><D2B5><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>,<2C><><EFBFBD>԰<EFBFBD><D4B0><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>м佻<D0BC><E4BDBB>
		double minX = resPoints[0].x;
		double maxX = resPoints[0].x;
		int minIndex = 0;
		int maxIndex = 0;

		for (int i=0;i<resPoints.size();++i) {
			if (resPoints[i].x < minX) minIndex = i;
			if (resPoints[i].x > maxX) maxIndex = i;
		}
		result.push_back(resPoints[minIndex]);
		result.push_back(resPoints[maxIndex]);
	}

	return resPoints.size();

}

// <20>ж<EFBFBD><D0B6><EFBFBD><EFBFBD><EFBFBD><EFBFBD>߶<EFBFBD><DFB6>Ƿ<EFBFBD><C7B7>ཻ
int TopologicalAnalysis::isLineIntersecting(const Line2D& l1, const Line2D& l2, Point2D& intersection) {
	float  numera = (l2.p2.x-l2.p1.x) * (l1.p1.y-l2.p1.y) - (l2.p2.y-l2.p1.y) * (l1.p1.x-l2.p1.x);
	float  numerb = (l1.p2.x-l1.p1.x) * (l1.p1.y-l2.p1.y) - (l1.p2.y-l1.p1.y) * (l1.p1.x-l2.p1.x);

	float  denom  = (l2.p2.y-l2.p1.y) * (l1.p2.x-l1.p1.x) - (l2.p2.x-l2.p1.x) * (l1.p2.y-l1.p1.y);

	if(denom == 0.0f)
	{
		if(numera == 0.0f && numerb == 0.0f)
		{
			intersection.x = l2.p2.x;
			intersection.y = l2.p2.y;
			return 2; //<2F>غ<EFBFBD>
		}
		return 3; //ƽ<><C6BD>
	}

	float ua = numera / denom;
	float ub = numerb / denom;

	if(ua >= 0.0f && ua <= 1.0f && ub >= 0.0f && ub <= 1.0f)
	{
		intersection.x = l1.p1.x + ua*(l1.p2.x - l1.p1.x);
		intersection.y = l1.p1.y + ua*(l1.p2.y - l1.p1.y);

/*
		double a1 = l1.p2.y - l1.p1.y;
		double b1 = l1.p1.x - l1.p2.x;
		double c1 = l1.p2.x * l1.p1.y - l1.p1.x * l1.p2.y;

		double a2 = l2.p2.y - l2.p1.y;
		double b2 = l2.p1.x - l2.p2.x;
		double c2 = l2.p2.x * l2.p1.y - l2.p1.x * l2.p2.y;
		double det = a1 * b2 - a2 * b1;

		intersection.x = (b2 * c1 - b1 * c2) / det;
		intersection.y = (a1 * c2 - a2 * c1) / det;*/

		return 1; //<2F>ཻ
	}

	return 0; //<2F><><EFBFBD>ཻ
}

//<2F><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С<EFBFBD><D0A1><EFBFBD><EFBFBD>
int TopologicalAnalysis::getLineWithPolygonMinWidth(const vector<double>& polygonX,const vector<double>& polygonY)
{
	int n = polygonX.size();
	double minWidth = INT64_MAX;
	int lineID = 0;

	for (int i=0;i<n-1;++i) //<2F><><EFBFBD><EFBFBD>ÿһ<C3BF><D2BB><EFBFBD><EFBFBD>
	{
		int startID = i;
		int endID = i + 1;
		double width = 0;
		for (int j=0;j<n-1;++j) //<2F><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>㵽ֱ<E3B5BD>ߵľ<DFB5><C4BE><EFBFBD>
		{
			if(startID == j || endID == j || (endID== n-1 && j==0)) continue;
			//<2F><><EFBFBD>㴹<EFBFBD><E3B4B9><EFBFBD><EFBFBD>
			double linePt1[2] = {polygonX[startID],polygonY[startID]};
			double linePt2[2] = {polygonX[endID],polygonY[endID]};
			double targetPt[2] = {0.0,0.0};
			double pt[2] = {polygonX[j],polygonY[j]};
			GetPointToLineVerticalCross(linePt1,linePt2,pt,targetPt);
			//<2F><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
			double dist;
			CalculateTwoPtsDistance(dist,pt[0],pt[1],targetPt[0],targetPt[1],3);
			if(dist > width) width = dist;
		}
		if(width < minWidth)
		{
			minWidth = width;
			lineID = startID;
		}
	}
	
	return lineID;
}

//<2F>жϵ<D0B6><CFB5><EFBFBD>ֱ<EFBFBD>ߵ<EFBFBD>λ<EFBFBD>ù<EFBFBD>ϵ
int TopologicalAnalysis::pointPosition(const Point2D& pt,const Line2D& line)
{
	double cross = crossProduct(line.p1, line.p2, pt);
	if (cross > 0) {
		return 1; //<2F><><EFBFBD><EFBFBD>
	} else if (cross < 0) {
		return -1; //<2F>Ҳ<EFBFBD>
	} else {
		return 0; //<2F><><EFBFBD><EFBFBD>
	}
}


// <20>жϵ<D0B6>P<EFBFBD>Ƿ<EFBFBD><C7B7>ڶ<EFBFBD><DAB6><EFBFBD><EFBFBD><EFBFBD>polygon<6F><6E>
bool TopologicalAnalysis::isPointInPolygon(const Point2D& P, vector<double>& regionLons,vector<double>& regionLats) {
	int n =  regionLons.size();
	bool inside = false;
	double xints;
	for(int i=0, j=n-1; i<n; j=i++) {
		double xi =  regionLons[i], yi = regionLats[i];
		double xj =  regionLons[j], yj = regionLats[j];

		// <20><><EFBFBD>ڶ<EFBFBD><DAB6><EFBFBD><EFBFBD>εĶ<CEB5><C4B6><EFBFBD><EFBFBD><EFBFBD>Ҳ<EFBFBD><D2B2>һ<EFBFBD><D2BB><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>򻯴<EFBFBD><F2BBAFB4><EFBFBD>δ<EFBFBD><CEB4><EFBFBD><EFBFBD>
		// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ƿ<EFBFBD><C7B7>뵱ǰ<EBB5B1><C7B0><EFBFBD>ཻ
		if(((yi > P.y) != (yj > P.y)) && // <20>ߵ<EFBFBD><DFB5><EFBFBD><EFBFBD>˵<EFBFBD><CBB5><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
			(P.x < (xj - xi) * (P.y - yi) / (yj - yi) + xi)) { // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ཻ<EFBFBD><E0BDBB><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
				inside = !inside; // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>һ<EFBFBD><D2BB><EFBFBD>ı<EFBFBD>inside״̬
		}
	}
	return inside;
}

//ֱ<><D6B1><EFBFBD><EFBFBD><EFBFBD>߶εĽ<CEB5><C4BD><EFBFBD>
bool TopologicalAnalysis::findIntersection(double x1, double y1, double x2, double y2, // ֱ<><D6B1>1<EFBFBD><31><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
	double x3, double y3, double x4, double y4,Point2D& intersection) { // <20>߶ε<DFB6><CEB5><EFBFBD><EFBFBD><EFBFBD><EFBFBD>˵<EFBFBD>
		// <20><><EFBFBD><EFBFBD>ֱ<EFBFBD><D6B1>1<EFBFBD>ķ<EFBFBD><C4B7><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
		double dx1 = x2 - x1;
		double dy1 = y2 - y1;

		// <20><><EFBFBD><EFBFBD><EFBFBD>߶εķ<CEB5><C4B7><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
		double dx2 = x4 - x3;
		double dy2 = y4 - y3;

		// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ж<EFBFBD><D0B6>Ƿ<EFBFBD>ƽ<EFBFBD>У<EFBFBD><D0A3><EFBFBD><EFBFBD>޽<EFBFBD><DEBD>㣩
		double cross = -dx1 * dy2 + dy1 * dx2;
		if (cross == 0) return false;

		// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>t1<74><31>ֱ<EFBFBD><D6B1>1<EFBFBD>ϵĵ㣩<C4B5><E3A3A9>t2<74><32><EFBFBD>߶<EFBFBD><DFB6>ϵĵ㣩
		double t1 = (dx2 * (y3 - y1) - dy2 * (x3 - x1)) / cross;
		double t2 = (dx1 * (y3 - y1) - dy1 * (x3 - x1)) / cross;

		// <20>ж<EFBFBD>t2<74>Ƿ<EFBFBD><C7B7><EFBFBD><EFBFBD>߶<EFBFBD><DFB6><EFBFBD>
		if (t2 >= 0 && t2 <= 1) {
			// <20><><EFBFBD>㲢<EFBFBD><E3B2A2><EFBFBD>ؽ<EFBFBD><D8BD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
			intersection.x = x1 + t1 * dx1;
			intersection.y = y1 + t1 * dy1;
			return true;
		} else {
			// <20><><EFBFBD><EFBFBD>t2<74><32><EFBFBD><EFBFBD>[0, 1]<5D><><EFBFBD><EFBFBD><EFBFBD>ڣ<EFBFBD>˵<EFBFBD><CBB5><EFBFBD><EFBFBD><EFBFBD>ཻ<EFBFBD><E0BDBB><EFBFBD>߶<EFBFBD><DFB6><EFBFBD>
			return false;
		}
}