分享一下之前做项目写的空间计算几何源码，
第一部分，涵盖简单计算几何方法的实现。例如：外包框计算、距离计算、面积与长度计算、空间关系的判断、方位角、交点计算、多边形中心点、多边形最大内接圆、凸包、投影点等。
第二部分主要是复杂的空间计算方法：多边形的交差并补计算、缓冲区、求平行线、求延长线等
第三部分主要是分享一下空间拓扑。
第二、三部分以后补上。
一. 头文件
//************************************** 1.距离计算 **************************************// //距离：点point到Geometry的最短距离（bsqrt表示是否开方） static double PtDistToGeometry(const Point& point, const Geometry* geometry, bool bsqrt = true); static double PtDistToGeometryP(const Point& point, const Geometry* geometry, Point& pointRet, bool bsqrt = true); //距离：两点之间的距离（bsqrt表示是否开方） static double PtDist(double x1,double y1,double x2,double y2, bool bsqrt = true); static double PtDist(const Point& point1, const Point& point2, bool bsqrt = true); //距离：点point到线段的最短距离（bsqrt表示是否开方）,并返回距该点最近的点 static double PtDistToSegment(double x, double y, double x1, double y1, double x2, double y2,double &XResult,double &YResult, bool bsqrt = true); static double PtDistToSegment(const Point& point, const Point& point1, const Point& point2, Point& pointRet, bool bsqrt = true); static double PtDistToSegment(double x, double y, double x1, double y1, double x2, double y2, bool bsqrt = true); static double PtDistToSegment(const Point& point, const Point& point1, const Point& point2, bool bsqrt = true); //距离：点point到折线LineString的最短距离（bsqrt表示是否开方） static double PtDistToLineString(const Point& point, const LineString* lineString, bool bsqrt = true); static double PtDistToLineString(const Point& point, const PointVec& pointVec, bool bsqrt = true); static double PtDistToLineStringP(const Point& point, const LineString* lineString, Point& pointRet, bool bsqrt = true);//返回最近点坐标 static double PtDistToLineStringP(const Point& point, const PointVec& pointVec, Point& pointRet, bool bsqrt = true);//返回最近点坐标 static double PtDistToLineStringP_SegIndex(const Point& point, const LineString* lineString, Point& pointRet, int& segIndex, bool bsqrt = true);//返回最近点坐标与结果点在折线上的序号 static double PtDistToLineStringP_SegIndex(const Point& point, const PointVec& pointVec, Point& pointRet, int& segIndex, bool bsqrt = true);//返回最近点坐标与结果点在折线上的序号 //距离：点point到多线MultiLineString的最短距离（bsqrt表示是否开方） static double PtDistToMultiLineString(const Point& point, const MultiLineString* multiLineString, bool bsqrt = true); static double PtDistToMultiLineStringP(const Point& point, const MultiLineString* multiLineString, Point& pointRet, bool bsqrt = true);//返回最近点坐标 //距离：点point到环LinearRing的最短距离（bsqrt表示是否开方, 在LinearRing内距离为负数） static double PtDistToLinearRing(const Point& point, const LinearRing* linearRing, bool bsqrt = true); static double PtDistToLinearRingP(const Point& point, const LinearRing* linearRing, Point& pointRet, bool bsqrt = true);//返回最近点坐标 static double PtDistToLinearRingP_SegIndex(const Point& point, const LinearRing* linearRing, Point& pointRet, int& segIndex, bool bsqrt = true);//返回最近点坐标与结果点在折线上的序号 //距离：点point到多边形Polygon边线的最短距离（bsqrt表示是否开方, 在Polygon内距离为负数） static double PtDistToPolygon(const Point& point, const Polygon* polygon, bool bsqrt = true); static double PtDistToPolygonP(const Point& point, const Polygon* polygon, Point& pointRet, bool bsqrt = true);//返回最近点坐标 static double PtDistToPolygonP_RingAndSegIndex(const Point& point, const Polygon* polygon, Point& pointRet, int& ringIndex, int& segIndex, bool bsqrt = true);//返回最近点坐标,环序号以及线段序号 //距离：点point到多面MultiPolygon边线的最短距离（bsqrt表示是否开方, 在MultiPolygon内距离为负数） static double PtDistToMultiPolygon(const Point& point, const MultiPolygon* multiPolygon, bool bsqrt = true); static double PtDistToMultiPolygonP(const Point& point, const MultiPolygon* multiPolygon, Point& pointRet, bool bsqrt = true);//返回最近点坐标 //距离：点到矩形的最短距离（bsqrt表示是否开方, 在矩形内距离为负数） static double PtDistToRectangle(const Point& point, const Envelope& envelope, bool bsqrt = true); static double PtDistToRectangleP(const Point& point, const Envelope& envelope, Point& pointRet, bool bsqrt = true);//返回最近点坐标 //距离：点到圆的最短距离（bsqrt表示是否开方, 在圆内距离为负数） static double PtDistToCircle(const Point& point, const Point& pointCenter, double radias, bool bsqrt = true); static double PtDistToCircleP(const Point& point, const Point& pointCenter, double radias, Point& pointRet, bool bsqrt = true);//返回最近点坐标 //距离：点到椭圆的最短距离（bsqrt表示是否开方, 在椭圆内距离为负数） static double PtDistToEllipse(const Point& point, const Point& pointCenter, double radiasA, double radiasB, bool bsqrt = true); static double PtDistToEllipseP(const Point& point, const Point& pointCenter, double radiasA, double radiasB, Point& pointRet, bool bsqrt = true);//返回最近点坐标 //***************************** 2.关系判断 ********************************// //关系判断(点2点)：判断两个点是否相等 static bool IsPointsEqual(const Point& point1,const Point& point2); //关系判断(点2线)：判断点是否在线段上 算法：(p在线段l所在的直线上) && (点p在以线段l为对角线的矩形内) static bool IsPointOnLine(const Point& point, const Point& point1, const Point& point2); static bool IsPointOnLine(const Point& point, const PointVec& pointVec); //判断点point是否在折线上 static bool IsPointOnLine(const Point& point, const LineString* lineString);//判断点point是否在折线lineString上 static bool IsPointOnLine(const Point& point, const PointVec& pointVec, int& segIndex);//判断点point是否在折线上，并返回该点在折线上的序号 static bool IsPointOnLine(const Point& point, const LineString* lineString, int& segIndex);//判断点point是否在折线lineString上，并返回该点在折线上的序号 //关系判断(点2线)：判断三点是否共线 static bool IsPointsInOneLine(const Point& point1, const Point& point2, const Point& point3); //关系判断(点2线)：点与线段的几何关系: -1左侧 0上 1右侧 2正向延长线上 3反向延长线上 static int PointInLine(const Point& point, const Point& point1, const Point& point2); //关系判断(点2面)：判断点与矩形、圆、椭圆、多边形之间的关系： -1 外 0 上 1 内 static int PointInRectangle(const Point& point, const Envelope& envelope); static int PointInCircle(const Point& point, const Point& pointCenter, double radias); static int PointInEllipse(const Point& point, const Point& pointCenter, double radiasA, double radiasB); static int PointInLinearRing(const Point& point, const LinearRing* linearRing); static int PointInPolygon(const Point& point, const Polygon* polygon); //关系判断(线2线)：线段u和线段v是否平行 static bool IsLinesParallel(const Point& u1, const Point& u2, const Point& v1, const Point& v2); //关系判断(线2线)：线段u和线段v是否相交(包括端点和部分重合) static bool IsLinesIntersect(const Point& u1, const Point& u2, const Point& v1, const Point& v2); static bool IsLinesIntersect(const Point& u1, const Point& u2, const PointVec& pointVec);//线段u和折线段是否相交 static bool IsLinesIntersect(const Point& u1, const Point& u2, const LineString* lineString);//线段u和折线段lineString是否相交 static bool IsLinesIntersect(const LineString* lineStringU, const LineString* lineStringV);//折线段lineStringU和折线段lineString是否相交 //关系判断(线2线)：线段u和线段v是否相交(不包括端点和部分重合) static bool IsLinesIntersect_Ignore_TwoSides(const Point& u1, const Point& u2, const Point& v1, const Point& v2); static bool IsLinesIntersect_Ignore_TwoSides(const Point& u1, const Point& u2, const PointVec& pointVec);//线段u和折线段是否相交 static bool IsLinesIntersect_Ignore_TwoSides(const Point& u1, const Point& u2, const LineString* lineString);//线段u和折线段lineString是否相交 static bool IsLinesIntersect_Ignore_TwoSides(const LineString* lineStringU, const LineString* lineStringV);//折线段lineStringU和折线段lineString是否相交 //关系判断(线2线)：线段u与直线v是否相交 static bool IsLinesIntersect_Seg_Line(const Point& u1, const Point& u2, const Point& v1, const Point& v2); //关系判断(线2线)：直线u与直线v是否相交 static bool IsLinesIntersect_Line_Line(const Point& u1, const Point& u2, const Point& v1, const Point& v2); //关系判断(线2面)：判断线段与矩形、圆、椭圆、多边形是否相交 static bool IsLineIntersectRectangle(const Point& point1, const Point& point2, const Envelope& envelope); static bool IsLineIntersectCircle(const Point& point1, const Point& point2, const Point& pointCenter, double radias); static bool IsLineIntersectEllipse(const Point& point1, const Point& point2, const Point& pointCenter, double radiasA, double radiasB); static bool IsLineIntersectPolygon(const Point& point1, const Point& point2, const Polygon* polygon); static bool IsLineIntersectRectangle(const LineString* lineString, const Envelope& envelope); static bool IsLineIntersectCircle(const LineString* lineString, const Point& pointCenter, double radias); static bool IsLineIntersectEllipse(const LineString* lineString, const Point& pointCenter, double radiasA, double radiasB); static bool IsLineIntersectPolygon(const LineString* lineString, const Polygon* polygon); //关系判断(线2面)：判断线段与矩形、圆、椭圆、多边形是否包含 static bool IsLineInRectangle(const Point& point1, const Point& point2, const Envelope& envelope); static bool IsLineInCircle(const Point& point1, const Point& point2, const Point& pointCenter, double radias); static bool IsLineInEllipse(const Point& point1, const Point& point2, const Point& pointCenter, double radiasA, double radiasB); static bool IsLineInPolygon(const Point& point1, const Point& point2, const Polygon* polygon); static bool IsLineInRectangle(const LineString* lineString, const Envelope& envelope); static bool IsLineInCircle(const LineString* lineString, const Point& pointCenter, double radias); static bool IsLineInEllipse(const LineString* lineString, const Point& pointCenter, double radiasA, double radiasB); static bool IsLineInPolygon(const LineString* lineString, const Polygon* polygon); //关系判断(面2面)：判断面与面是否相交 static bool IsRectangleIntersectRectangle(const Envelope& envelope1, const Envelope& envelope2); static bool IsRectangleIntersectCircle(const Envelope& envelope, const Point& pointCenter, double radias); static bool IsRectangleIntersectPolygon(const Envelope& envelope, const Polygon* polygon); static bool IsCircleIntersectCircle(const Point& pointCenter1, double radias1, const Point& pointCenter2, double radias2); static bool IsCircleIntersectPolygon(const Point& pointCenter, double radias, const Polygon* polygon); static bool IsPolygonIntersectPolygon(const Polygon* polygon1, const Polygon* polygon2); //关系判断(面2面)：判断面与面是否包含 static bool IsRectangleInRectangle(const Envelope& envelope1, const Envelope& envelope2); static bool IsRectangleInCircle(const Envelope& envelope, const Point& pointCenter, double radias); static bool IsRectangleInPolygon(const Envelope& envelope, const Polygon* polygon); static bool IsCircleInRectangle(const Point& pointCenter, double radias, const Envelope& envelope); static bool IsCircleInCircle(const Point& pointCenter1, double radias1, const Point& pointCenter2, double radias2); static bool IsCircleInPolygon(const Point& pointCenter, double radias, const Polygon* polygon); static bool IsPolygonInRectangle(const Polygon* polygon, const Envelope& envelope); static bool IsPolygonInCircle(const Polygon* polygon, const Point& pointCenter, double radias); static bool IsPolygonInPolygon(const Polygon* polygon1, const Polygon* polygon2); //***************************** 3.求值 ********************************// //投影点：计算点P0(x，y)在线上的投影点（垂点）,并返回距投影点在线段上的比例：在正向延长线上值应大于1， 在线段上为0-1， 在反向延长线上<0 static double PointPrj2Line(double x, double y, double x1, double y1, double x2, double y2, double &XResult, double &YResult); static double PointPrj2Line(const Point& point, const Point& point1, const Point& point2, Point& pointRet); //旋转点：（pointLoc为定位点, pointOld为旋转前坐标，xNew与yNew为旋转后坐标，angle为旋转角度） static void PointRotate(double xLoc, double yLoc, double xOld, double yOld, double& xNew, double& yNew, double angle); static Point PointRotate(const Point& pointLoc, const Point& pointOld, double angle); //分割点：计算线上按比例分割的点 static void PointOnLine(double x1, double y1, double x2, double y2, double dScale, double &XResult,double &YResult); static Point PointOnLine(const Point& point1, const Point& point2, double dScale); static Point PointOnLine(const LineString* lineString, double dScale); static Point PointOnLine(const PointVec& pointVec, double dScale); static double ScaleOnLine(const Point& point, const PointVec& pointVec);//返回点point在折线上的比例 static double ScaleOnLine(const Point& point, const LineString* lineString);//返回点point在折线lineString上的比例 //中点：计算线段中点 static void PointCenter(double x1, double y1, double x2, double y2, double& x, double& y); static Point PointCenter(const Point& point1, const Point& point2); //延伸点：以point2为起始延伸点根据延长距离计算延伸点（extendDistance ：延伸距离，负数表示反方向延伸） static void PointExtend(double x1, double y1, double x2, double y2, double extendDistance, double &XResult,double &YResult); static Point PointExtend(const Point& point1, const Point& point2, double extendDistance); static Point PointExtend(const LineString* lineString, double extendDistance, bool beginStart = false); //beginStart:从第一个点开始延伸 static Point PointExtend(const PointVec& pointVec, double extendDistance, bool beginStart = false); //垂直点：过直线（point1， point2）上一点point，求垂直于该点的直线外一点(直线左手方向distance为正数) static void PPLDN(double x,double y,double x1,double y1,double x2, double y2,double distance, double& resultx,double& resulty); static void PPLDN(const Point& point, const Point& point1, const Point& point2, Point& pointResult, double distance); //对称点：已知两点（point1， point2）的直线和第三点point3，求第三点的对称点pointResult static void P2LSymmetry(const Point& point1, const Point& point2, const Point& point3, Point& pointResult); static void P2LSymmetry(double x1,double y1,double x2,double y2, double x3, double y3, double &XResult,double &YResult); //对称点：求点point关于点symmetryPoint的对称点 static void P2PSymmetry(double x,double y,double xRef,double yRef, double& resultx,double& resulty); static void P2PSymmetry(const Point& point, const Point& symmetryPoint, Point& pointResult); //求交点：(如果两线平行且重合，取重合部分的两个端点) static bool GetIntersectPoint(const Point& u1, const Point& u2, const Point& v1, const Point& v2, PointVec& pointVecIntersect); static bool GetIntersectPoint(const Point& u1, const Point& u2, const LineString* lineString, PointVec& pointVecIntersect); //求交点：(如果两线平行且重合，交点是u1和u2) static bool GetIntersectPoint_Seg_Line(const Point& u1, const Point& u2, const Point& v1, const Point& v2, PointVec& pointVecIntersect); //求交点：(如果两线平行且重合，没有交点) static bool GetIntersectPoint_Line_Line(const Point& u1, const Point& u2, const Point& v1, const Point& v2, Point& pointRet); //求线段与矩形、交点 static void GetIntersectPoint_Seg_Rectangle(const Point& point1, const Point& point2, const Envelope& envelope, PointVec& pointVecIntersect); static void GetIntersectPoint_Seg_Circle(const Point& point1, const Point& point2, const Point& pointCenter, double radias, PointVec& pointVecIntersect); static void GetIntersectPoint_Seg_Ellipse(const Point& point1, const Point& point2, const Point& pointCenter, double radiasA, double radiasB, PointVec& pointVecIntersect); static void GetIntersectPoint_Seg_Polygon(const Point& point1, const Point& point2, const Polygon* polygon, PointVec& pointVecIntersect); //返回交点与交点所在的环序号以及环内第几段 //linearingIndexAndSegIndex：使用Point保存序号，X坐标表示环的序号，-1表示最外环；Y坐标表示当前环上的第几段，从0开始 static void GetIntersectPoint_Seg_Polygon(const Point& point1, const Point& point2, const Polygon* polygon, PointVec& pointVecIntersect, PointVec& linearingIndexAndSegIndex); //多边形或环上的某一点的凹凸性 -1凹 0共线 1凸 static int PointConvexity(Point pointCur, Point ptFirst, Point ptLast, bool isClockWise); static int PointConvexity(int pointIndex, const LinearRing* linearRing); static int PointConvexity(int pointIndex, const Polygon* polygon, int lineRingIndex = -1); //线上的某一点的凹凸性 -1凹 0共线 1凸：本来LineString没有凹凸之分，这里约定第三个点在前两个点组成向量的顺时针方向，则为凸；否则为凹 static int PointConvexity(Point pointCur, Point ptFirst, Point ptLast); static int PointConvexity(int pointIndex, const LineString* lineString); //计算点集的凸包 static void PointsConvexHull(const PointVec& points, PointVec& pointVecRet); static Polygon* PointsConvexHull(const PointVec& points); //计算点集的中心点 static Point PointsCentroid(const PointVec& points); //快速获取多边形内任意一点 static Point PolygonInnerPoint(const Polygon* polygon); //多边形重心 static Point LinearRingCentroid(const LinearRing* linearRing); static Point PolygonCentroid(const Polygon* polygon); //多边形最大内接圆（近似最大） static Point LinearRingMaxInnerCircle(const LinearRing* linearRing, double& radias, int iterateCount = 5); static Point PolygonMaxInnerCircle(const Polygon* polygon, double& radias, int iterateCount = 5); //方位角：线段的方位角（弧度0-2pi）, 起始方位：3点钟方向，逆时针 static double LineOrientation(double x1, double y1, double x2, double y2); static double LineOrientation(const Point& point1, const Point& point2); //夹角：以点point1、anglePoint、point2三点的夹角 static double LinesAngle(const Point& anglePoint, const Point& point1, const Point& point2); //夹角：前两点形成向量和后两点形成向量的夹角 static double LinesAngle(const Point& startPoint1, const Point& endPoint1, const Point& startPoint2, const Point& endPoint2); //角平分线：第一个向量顺时针或逆时针方向的角平分线 static Point LinesHalfAnglePoint(const Point& anglePoint, const Point& point1, const Point& point2, bool clockWise); //面积与长度 static double LineLength(const LineString* lineString, int ptFromIndex=0, int ptToIndex=ConstValue::MAXVALUE); static double LineLength(const PointVec& pointVec, int ptFromIndex=0, int ptToIndex=ConstValue::MAXVALUE); static double PolygonLength(const Polygon* polygon); static double MultiPolygonLength(const MultiPolygon* multiPolygon); static double MultiPolygonArea(const MultiPolygon* multiPolygon); static double TriangleAera(const Point& point1, const Point& point2, const Point& point3);//三角形面积 static double PointVecAera(const PointVec& pointVec);//点数组面积 //外包框 static Envelope MakeEnvelope(const Point& point1, const Point& point2, double unvalidBuffer = 0.0000001); static Envelope MakeEnvelope(const Point& point1, double buffer);
二.Cpp文件
/*r=multx(sp,ep,op),得到(sp-op)*(ep-op)的叉积 r>0:ep在矢量opsp的逆时针方向； r=0：opspep三点共线； r<0:ep在矢量opsp的顺时针方向*/ inline double multx(const Point& p1, const Point& p2, const Point& p0) { return (p1.X() - p0.X()) * (p2.Y() - p0.Y()) - (p2.X() - p0.X()) * (p1.Y() - p0.Y()); } /*r=multd(p1,p2,op),得到矢量(p1-op)和(p2-op)的点积，如果两个矢量都非零矢量 r<0:两矢量夹角为锐角；r=0：两矢量夹角为直角；r>0:两矢量夹角为钝角*/ inline double multd(const Point& p1, const Point& p2, const Point& p0) { return ((p1.X()-p0.X())*(p2.X()-p0.X()) + (p1.Y()-p0.Y())*(p2.Y()-p0.Y())); } inline bool getIntersectPoint(const Point& u1, const Point& u2, const Point& v1, const Point& v2, Point& pointRet) { double t = ((u1.X()-v1.X())*(v1.Y()-v2.Y())-(u1.Y()-v1.Y())*(v1.X()-v2.X())) / ((u1.X()-u2.X())*(v1.Y()-v2.Y())-(u1.Y()-u2.Y())*(v1.X()-v2.X())); pointRet.setX(u1.X() + (u2.X() - u1.X()) * t); pointRet.setY(u1.Y() + (u2.Y() - u1.Y()) * t); return true; } inline double ptDist(double x1,double y1,double x2,double y2, bool bsqrt = true) { double dis = (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1); return bsqrt ? sqrt(dis) : dis ; } inline double ptDist(const Point& point1, const Point& point2, bool bsqrt = true) { return ptDist(point1.getX(), point1.getY(), point2.getX(), point2.getY(), bsqrt); } inline void ptCenter(double x1,double y1,double x2,double y2, double& x, double& y) { x = (x2 + x1) * 0.5; y = (y2 + y1) * 0.5; } inline void uniquePoints(std::vector & points) { std::sort(points.begin(), points.end(),[&](const Point& l, const Point& r) { return l.X() == r.X() ? l.Y() < r.Y() : l.X() < r.X(); });//对所有点进行排序 按x坐标升序排序， 当x相等时取y小的 points.erase(unique(points.begin(), points.end()), points.end()); } inline void uniquePoints(GeoVectorT & points) { std::vector points2; points2.assign(points.begin(), points.end()); uniquePoints(points2); points.clear(); for(size_t i = 0; i < points2.size(); i++) points.push_back(points2[i]); } inline Point getNearistPoint(Point point, GeoVectorT & points, bool excludeSelf = false) { Point ptRet(ConstValue::MAXVALUE, ConstValue::MAXVALUE); double dmin = ConstValue::MAXVALUE; for(size_t j = 0; j < points.size(); j++) { if(excludeSelf && GeoSpatialCaculation::IsPointsEqual(point, points[j])) continue; double d = GeoSpatialCaculation::PtDist(point, points[j], false); if(dmin > d) { dmin = d; ptRet = points[j]; } } return ptRet; } inline void lineIntersectEllipseP(const Point& point1, const Point& point2, const Point& pointCenter, double radiasA, double radiasB, GeoVectorT & pointVecIntersect) { //利用直线方程斜截式计算椭圆与直线的交点 double x1 = point1.getX(); double y1 = point1.getY(); double x2 = point2.getX(); double y2 = point2.getY(); double rx1 = ConstValue::MAXVALUE; double rx2 = ConstValue::MAXVALUE; double ry1 = ConstValue::MAXVALUE; double ry2 = ConstValue::MAXVALUE; //直线方程表达式为x = C(C为常数); if(x1 == x2) { rx1 = x1; rx2 = x1; ry1 = sqrt(radiasB * radiasB * (radiasA * radiasA - x1 * x1) / (radiasA * radiasA)); ry2 = -1 * ry1; } else if(y1 == y2)//直线方程表达式为y = C(C为常数); { ry1 = y1; ry2 = y1; rx1 = sqrt(radiasA * radiasA * (radiasB * radiasB - y1 * y1) / (radiasB * radiasB)); rx2 = -1 * rx1; } else//直线方程表达式为y = k*x + C(C为常数); { double k = (y2 - y1)/(x2 - x1); double c = y2 - k * x2; //判别式b*b-4*a*c double tmp = sqrt((4 * k * k * c * c * radiasA * radiasA * radiasA * radiasA) - 4 * radiasA * radiasA * (radiasB * radiasB + radiasA * radiasA * k * k) * (c * c - radiasB * radiasB)); rx1 = (-2 * k * c * radiasA * radiasA + tmp) / (2 * (radiasB * radiasB + radiasA * radiasA * k * k)); rx2 = (-2 * k * c * radiasA * radiasA - tmp) / (2 * (radiasB * radiasB + radiasA * radiasA * k * k)); ry1 = k * rx1 + c; ry2 = k * rx2 + c; } Point rp1 = Point(rx1, ry1); Point rp2 = Point(rx2, ry2); bool b1 = GeoSpatialCaculation::IsPointOnLine(rp1, point1, point2); bool b2 = GeoSpatialCaculation::IsPointOnLine(rp2, point1, point2); if(b1) pointVecIntersect.push_back(rp1); if(b2) pointVecIntersect.push_back(rp2); } inline PointVector getCentralPoints(double xMin, double yMin, double xMax, double yMax, int iterateCount, bool full = true) { PointVector pointVector; for (int i = 0; i < iterateCount; i++) { if(full == false && i != iterateCount - 1) continue; int m = pow(2.0, i); int n = pow(2.0, (i + 1)); for (int j = m - 1; j >= 0; j--) { double y = yMax + (yMin - yMax)/n*(2*j + 1); for (int k = 0; k < m; k++) { double x = xMin + (xMax - xMin)/n*(2*k + 1); Point p = Point(x, y); pointVector.push_back(p); } } } return pointVector; } inline PointVector getEnvelopePoints(Envelope envelope) { PointVector pts; pts.push_back(Point(envelope.getMinX(), envelope.getMaxY())); pts.push_back(Point(envelope.getMaxX(), envelope.getMaxY())); pts.push_back(Point(envelope.getMaxX(), envelope.getMinY())); pts.push_back(Point(envelope.getMinX(), envelope.getMinY())); pts.push_back(Point(envelope.getMinX(), envelope.getMaxY())); return pts; } void getIntersectPoint_Seg_Polygon(const Point& point1, const Point& point2, const Polygon* polygon, GeoVectorT & pointVecIntersect, GeoVectorT & linearingIndexAndSegIndex, bool returnIndexVec = false) { pointVecIntersect.clear(); if(returnIndexVec) linearingIndexAndSegIndex.clear(); Envelope ev = GeoSpatialCaculation::MakeEnvelope(point1, point2); if(!ev.intersect(polygon->getEnvelope())) return; int n = (int)polygon->getInteriorRingsCount(); for(int i = -1; i < n; i++) { const LinearRing* line = (i == -1) ? polygon->getExteriorRing() : polygon->getInteriorRing(i); if(line == 0) continue; if(!ev.intersect(line->getEnvelope())) continue; for(size_t j = 1; j < line->getPointsCount(); j++) { Point pt0 = line->getPoint(j-1); Point pt1 = line->getPoint(j); //printf("%d \n", j); GeoVectorT pts; bool bIntersect = GeoSpatialCaculation::GetIntersectPoint(point1, point2, pt0, pt1, pts); if(!bIntersect) continue; pointVecIntersect.insert(pointVecIntersect.end(), pts.begin(), pts.end()); if(returnIndexVec) { for(size_t k = 0; k < pts.size(); k++) linearingIndexAndSegIndex.push_back(Point(i, j-1)); } } } //uniquePoints(pointVecIntersect); }
//************************************** 1.距离计算 **************************************//double GeoSpatialCaculation::PtDistToGeometry(const Point& point, const Geometry* geometry, bool bsqrt){ if(geometry == 0) return ConstValue::MAXVALUE; GeometryType geometryType = geometry->getGeometryType(); switch(geometryType) { case GeometryType_Point: return PtDist(point.getX(), point.getY(), ((const Point*)geometry)->getX(), ((const Point*)geometry)->getY(), bsqrt); case GeometryType_LineString: return PtDistToLineString(point, (const LineString*)geometry, bsqrt); case GeometryType_LinearRing: return PtDistToLinearRing(point, (const LinearRing*)geometry, bsqrt); case GeometryType_Polygon: return PtDistToPolygon(point, (const Polygon*)geometry, bsqrt); case GeometryType_MultiLineString: return PtDistToMultiLineString(point, (const MultiLineString*)geometry, bsqrt); case GeometryType_MultiPolygon: return PtDistToMultiPolygon(point, (const MultiPolygon*)geometry, bsqrt); case GeometryType_GeometryCollection: { const GeometryCollection* gc = (const GeometryCollection*)geometry; double minDis = ConstValue::MAXVALUE; for (size_t j = 0; j < gc->getCount(); j++) { double dis = PtDistToGeometry(point, gc->getGeometry(j), bsqrt); if(minDis > dis) minDis = dis; } return minDis; } } return ConstValue::MAXVALUE;}double GeoSpatialCaculation::PtDistToGeometryP(const Point& point, const Geometry* geometry, Point& pointRet, bool bsqrt){ if(geometry == 0) return ConstValue::MAXVALUE; GeometryType geometryType = geometry->getGeometryType(); switch(geometryType) { case GeometryType_Point: { pointRet = point; return PtDist(point.getX(), point.getY(), ((const Point*)geometry)->getX(), ((const Point*)geometry)->getY(), bsqrt); } case GeometryType_LineString: return PtDistToLineStringP(point, (const LineString*)geometry, pointRet, bsqrt); case GeometryType_LinearRing: return PtDistToLinearRingP(point, (const LinearRing*)geometry, pointRet, bsqrt); case GeometryType_Polygon: return PtDistToPolygonP(point, (const Polygon*)geometry, pointRet, bsqrt); case GeometryType_MultiLineString: return PtDistToMultiLineStringP(point, (const MultiLineString*)geometry, pointRet, bsqrt); case GeometryType_MultiPolygon: return PtDistToMultiPolygonP(point, (const MultiPolygon*)geometry, pointRet, bsqrt); case GeometryType_GeometryCollection: { const GeometryCollection* gc = (const GeometryCollection*)geometry; double minDis = ConstValue::MAXVALUE; for (size_t j = 0; j < gc->getCount(); j++) { double dis = PtDistToGeometryP(point, gc->getGeometry(j), pointRet, bsqrt); if(minDis > dis) minDis = dis; } return minDis; } } return ConstValue::MAXVALUE;}double GeoSpatialCaculation:: PtDist(double x1,double y1,double x2,double y2, bool bsqrt/* = true*/){ return ptDist(x1,y1,x2,y2, bsqrt);}double GeoSpatialCaculation:: PtDist(const Point& point1, const Point& point2, bool bsqrt/* = true*/){ return ptDist(point1, point2, bsqrt);}double GeoSpatialCaculation:: PtDistToSegment(double x, double y, double x1, double y1, double x2, double y2,double &XResult,double &YResult, bool bsqrt/* = true*/){ double cross = (x2 - x1) * (x - x1) + (y2 - y1) * (y - y1); if (cross <= 0) { XResult = x1; YResult = y1; return ptDist(x,y, XResult,YResult, bsqrt); } double d2 = (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1); if (cross >= d2) { XResult = x2; YResult = y2; return ptDist(x,y, XResult,YResult, bsqrt); } double r = cross / d2; XResult = x1 + (x2 - x1) * r; YResult = y1 + (y2 - y1) * r; return ptDist(x,y, XResult,YResult, bsqrt);}double GeoSpatialCaculation:: PtDistToSegment(const Point& point, const Point& point1, const Point& point2, Point& pointRet, bool bsqrt/* = true*/){ double x = 0; double y = 0; double dis = PtDistToSegment(point.X(), point.Y(), point1.X(), point1.Y(), point2.X(), point2.Y(), x, y, bsqrt); pointRet.setX(x); pointRet.setY(y); return dis;}double GeoSpatialCaculation:: PtDistToSegment(double x, double y, double x1, double y1, double x2, double y2, bool bsqrt/* = true*/){ double retx,rety = 0; return PtDistToSegment(x, y, x1, y1, x2, y2, x, y, bsqrt);}double GeoSpatialCaculation:: PtDistToSegment(const Point& point, const Point& point1, const Point& point2, bool bsqrt/* = true*/){ return PtDistToSegment(point.getX(), point.getY(), point1.getX(), point1.getY(), point2.getX(), point2.getY(), bsqrt);}double GeoSpatialCaculation:: PtDistToLineString(const Point& point, const LineString* lineString, bool bsqrt/* = true*/){ return PtDistToLineString(point, lineString->getPoints(), bsqrt);}double GeoSpatialCaculation:: PtDistToLineString(const Point& point, const PointVec& pointVec, bool bsqrt /*= true*/){ Point Pt; return PtDistToLineStringP(point, pointVec, Pt, bsqrt);}double GeoSpatialCaculation:: PtDistToLineStringP(const Point& point, const LineString* lineString, Point& pointRet, bool bsqrt/* = true*/){ return PtDistToLineStringP(point, lineString->getPoints(), pointRet, bsqrt);}double GeoSpatialCaculation:: PtDistToLineStringP(const Point& point, const PointVec& pointVec, Point& pointRet, bool bsqrt/* = true*/){ int segIndex = 0; return PtDistToLineStringP_SegIndex(point, pointVec, pointRet, segIndex, bsqrt);}double GeoSpatialCaculation:: PtDistToLineStringP_SegIndex(const Point& point, const LineString* lineString, Point& pointRet, int& segIndex, bool bsqrt/* = true*/){ return PtDistToLineStringP_SegIndex(point, lineString->getPoints(), pointRet, segIndex, bsqrt);}double GeoSpatialCaculation:: PtDistToLineStringP_SegIndex(const Point& point, const PointVec& pointVec, Point& pointRet, int& segIndex, bool bsqrt/* = true*/){ double minDis = ConstValue::MAXVALUE; Point Pt; for (size_t j = 1; j < pointVec.size(); j++) { double dis = PtDistToSegment(point, pointVec[j-1], pointVec[j], Pt, false); if(minDis > dis) { minDis = dis; pointRet = Pt; segIndex = j - 1; } if(dis == 0) break; } return bsqrt ? sqrt(minDis) : minDis;}double GeoSpatialCaculation:: PtDistToMultiLineString(const Point& point, const MultiLineString* multiLineString, bool bsqrt/* = true*/){ Point Pt; return PtDistToMultiLineStringP(point, multiLineString, Pt, bsqrt);}double GeoSpatialCaculation:: PtDistToMultiLineStringP(const Point& point, const MultiLineString* multiLineString, Point& pointRet, bool bsqrt/* = true*/){ double minDis = ConstValue::MAXVALUE; Point Pt; for (size_t j = 0; j < multiLineString->getCount(); j++) { const LineString* obj = multiLineString->getLineString(j); double dis = PtDistToLineStringP(point, obj, Pt, false); if(minDis > dis) { minDis = dis; pointRet = Pt; } if(dis == 0) break; } return bsqrt ? sqrt(minDis) : minDis;}double GeoSpatialCaculation:: PtDistToLinearRing(const Point& point, const LinearRing* linearRing, bool bsqrt/* = true*/){ Point Pt; return PtDistToLinearRingP(point, linearRing, Pt, bsqrt);}double GeoSpatialCaculation:: PtDistToLinearRingP(const Point& point, const LinearRing* linearRing, Point& pointRet, bool bsqrt/* = true*/){ int segIndex = -1; return PtDistToLinearRingP_SegIndex(point, linearRing, pointRet, segIndex, bsqrt);}double GeoSpatialCaculation:: PtDistToLinearRingP_SegIndex(const Point& point, const LinearRing* linearRing, Point& pointRet, int& segIndex, bool bsqrt/* = true*/){ double d = PtDistToLineStringP_SegIndex(point, linearRing, pointRet, segIndex, bsqrt); if(d == 0) return 0; return linearRing->containPoint(point) ? -d : d;}double GeoSpatialCaculation:: PtDistToPolygon(const Point& point, const Polygon* polygon, bool bsqrt/* = true*/){ Point Pt; return PtDistToPolygonP(point, polygon, Pt, bsqrt);}double GeoSpatialCaculation:: PtDistToPolygonP(const Point& point, const Polygon* polygon, Point& pointRet, bool bsqrt/* = true*/){ int ringIndex,segIndex = 0; return PtDistToPolygonP_RingAndSegIndex(point, polygon, pointRet, ringIndex, segIndex, bsqrt);}double GeoSpatialCaculation:: PtDistToPolygonP_RingAndSegIndex(const Point& point, const Polygon* polygon, Point& pointRet, int& ringIndex, int& segIndex, bool bsqrt){ Point Pt; bool bInPoly = false; double minDis = ConstValue::MAXVALUE; int n = (int)polygon->getInteriorRingsCount(); for (int j = -1; j < n; j++) { const LinearRing* linearRing = (j == -1) ? polygon->getExteriorRing() : polygon->getInteriorRing(j); int segIndexLine = 0; double dis = PtDistToLinearRingP_SegIndex(point, linearRing, Pt, segIndexLine, false); double dis2 = fabs(dis); if(minDis > dis2) { minDis = dis2; pointRet = Pt; if(j != -1) bInPoly = dis >= 0; else bInPoly = dis < 0; ringIndex = j; segIndex = segIndexLine; } if(dis == 0) break; } double d = bsqrt ? sqrt(minDis) : minDis; return bInPoly ? -d : d;}double GeoSpatialCaculation:: PtDistToMultiPolygon(const Point& point, const MultiPolygon* multiPolygon, bool bsqrt/* = true*/){ Point Pt; return PtDistToMultiPolygonP(point, multiPolygon, Pt, bsqrt);}double GeoSpatialCaculation:: PtDistToMultiPolygonP(const Point& point, const MultiPolygon* multiPolygon, Point& pointRet, bool bsqrt/* = true*/){ Point Pt; double minDis = ConstValue::MAXVALUE; bool bInPoly = false; for (size_t j = 0; j < multiPolygon->getCount(); j++) { const Polygon* obj = multiPolygon->getPolygon(j); double dis = PtDistToPolygonP(point, obj, Pt, false); double dis2 = fabs(dis); if(minDis > dis2) { minDis = dis2; pointRet = Pt; bInPoly = dis < 0; } if(dis == 0) break; } minDis = bsqrt ? sqrt(minDis) : minDis; return bInPoly ? -minDis : minDis;}double GeoSpatialCaculation:: PtDistToRectangle(const Point& point, const Envelope& envelope, bool bsqrt /*= true*/){ Point pointRet; return PtDistToRectangleP(point, envelope, pointRet, bsqrt);}double GeoSpatialCaculation:: PtDistToRectangleP(const Point& point, const Envelope& envelope, Point& pointRet, bool bsqrt /*= true*/){ double d = PtDistToLineStringP(point, getEnvelopePoints(envelope), pointRet, bsqrt); if(d == 0) return 0; return PointInRectangle(point, envelope) ? -d : d;}double GeoSpatialCaculation:: PtDistToCircle(const Point& point, const Point& pointCenter, double radias, bool bsqrt /*= true*/){ Point pointRet; return PtDistToCircleP(point, pointCenter, radias, pointRet, bsqrt);}double GeoSpatialCaculation:: PtDistToCircleP(const Point& point, const Point& pointCenter, double radias, Point& pointRet, bool bsqrt /*= true*/){ if(point == pointCenter) { pointRet = Point(pointCenter.getX() - radias, pointCenter.getY()); return bsqrt ? -radias : -radias*radias; } double angle = LineOrientation(pointCenter, point); pointRet = Point(pointCenter.getX() + radias*cos(angle), pointCenter.getY() + radias*sin(angle)); return PointInCircle(point, pointCenter, radias) > 0 ? -ptDist(pointRet, point, bsqrt) : ptDist(pointRet, point, bsqrt);}double GeoSpatialCaculation:: PtDistToEllipse(const Point& point, const Point& pointCenter, double radiasA, double radiasB, bool bsqrt /*= true*/){ Point pointRet; return PtDistToEllipseP(point, pointCenter, radiasA, radiasB, pointRet, bsqrt);}double GeoSpatialCaculation:: PtDistToEllipseP(const Point& point, const Point& pointCenter, double radiasA, double radiasB, Point& pointRet, bool bsqrt /*= true*/){ if(point == pointCenter) { if(radiasA <= radiasB) pointRet = Point(pointCenter.getX() + radiasA, pointCenter.getY()); else pointRet = Point(pointCenter.getX(), pointCenter.getY() + radiasB); double d = radiasA < radiasB ? radiasA : radiasB; return bsqrt ? -d : -d * d; } double angle = LineOrientation(pointCenter, point); double p = sqrt((radiasA * radiasA * radiasB * radiasB) / (cos(angle)*cos(angle) * radiasB * radiasB + sin(angle)*sin(angle) * radiasA * radiasA)); pointRet = Point(pointCenter.getX() + p*cos(angle), pointCenter.getY() + p*sin(angle)); return PointInEllipse(point, pointCenter, radiasA, radiasB) >= 0 ? -ptDist(pointRet, point, bsqrt) : ptDist(pointRet, point, bsqrt);}//***************************** 2.关系判断 ********************************//bool GeoSpatialCaculation:: IsPointOnLine(const Point& point, const Point& point1, const Point& point2){ if(fabs(multx(point2, point, point1)) > ConstValue::ZERO) return false; if((point1.X() - point.X()) * (point2.X() - point.X()) <= 0 && (point1.Y() - point.Y()) * (point2.Y() - point.Y()) <= 0) return true; return false;}bool GeoSpatialCaculation:: IsPointOnLine(const Point& point, const PointVec& pointVec){ int index = 0; return IsPointOnLine(point, pointVec, index);}bool GeoSpatialCaculation:: IsPointOnLine(const Point& point, const LineString* lineString){ return IsPointOnLine(point, lineString->getPoints());}bool GeoSpatialCaculation:: IsPointOnLine(const Point& point, const PointVec& pointVec, int& segIndex){ segIndex = 0; int n = (int)pointVec.size(); if(n <= 1) return false; for(int i = 1; i < n; i++) { if(IsPointOnLine(point, pointVec[i-1], pointVec[i])) return true; segIndex++; } return false;}bool GeoSpatialCaculation:: IsPointOnLine(const Point& point, const LineString* lineString, int& segIndex){ return IsPointOnLine(point, lineString->getPoints(), segIndex);}bool GeoSpatialCaculation:: IsPointsInOneLine(const Point& point1, const Point& point2, const Point& point3){ return fabs(multx(point1, point2, point3)) < ConstValue::ZERO;}bool GeoSpatialCaculation:: IsPointsEqual(const Point& point1,const Point& point2) { return (fabs(point1.X()-point2.X()) < ConstValue::ZERO) && (fabs(point1.Y()-point2.Y()) < ConstValue::ZERO);}bool GeoSpatialCaculation:: IsLinesParallel(const Point& u1, const Point& u2, const Point& v1, const Point& v2){ return fabs((u1.X()-u2.X())*(v1.Y()-v2.Y())-(v1.X()-v2.X())*(u1.Y()-u2.Y())) < ConstValue::ZERO;}bool GeoSpatialCaculation:: IsLinesIntersect(const Point& u1, const Point& u2, const Point& v1, const Point& v2){ return((max(u1.X(),u2.X()) >= min(v1.X(),v2.X()))&& //排斥实验 (max(v1.X(),v2.X()) >= min(u1.X(),u2.X()))&& (max(u1.Y(),u2.Y()) >= min(v1.Y(),v2.Y()))&& (max(v1.Y(),v2.Y()) >= min(u1.Y(),u2.Y()))&& (multx(v1,u2,u1) * multx(u2,v2,u1)>=0)&& //跨立实验 (multx(u1,v2,v1) * multx(v2,u2,v1)>=0));}bool GeoSpatialCaculation:: IsLinesIntersect(const Point& u1, const Point& u2, const PointVec& pointVec){ if(pointVec.size() <= 1) return false; for (size_t j = 1; j < pointVec.size(); j++) { if(IsLinesIntersect(u1, u2, pointVec[j-1], pointVec[j])) return true; } return false;}bool GeoSpatialCaculation:: IsLinesIntersect(const Point& u1, const Point& u2, const LineString* lineString){ return IsLinesIntersect(u1, u2, lineString->getPoints());}bool GeoSpatialCaculation:: IsLinesIntersect(const LineString* lineStringU, const LineString* lineStringV){ if(lineStringU->getPointsCount() <= 1) return false; for (size_t j = 1; j < lineStringU->getPointsCount(); j++) { if(IsLinesIntersect(lineStringU->getPoint(j-1), lineStringU->getPoint(j), lineStringV)) return true; } return false;}bool GeoSpatialCaculation:: IsLinesIntersect_Ignore_TwoSides(const Point& u1, const Point& u2, const Point& v1, const Point& v2){ return((IsLinesIntersect(u1, u2, v1, v2))&& (!IsPointOnLine(v1, u1, u2))&& (!IsPointOnLine(v2, u1, u2))&& (!IsPointOnLine(u2, v1, v2))&& (!IsPointOnLine(u1, v1, v2)));}bool GeoSpatialCaculation:: IsLinesIntersect_Ignore_TwoSides(const Point& u1, const Point& u2, const PointVec& pointVec){ if(pointVec.size() <= 1) return false; for (size_t j = 1; j < pointVec.size(); j++) { if(IsLinesIntersect_Ignore_TwoSides(u1, u2, pointVec[j-1], pointVec[j])) return true; } return false;}bool GeoSpatialCaculation:: IsLinesIntersect_Ignore_TwoSides(const Point& u1, const Point& u2, const LineString* lineString){ return IsLinesIntersect_Ignore_TwoSides(u1, u2, lineString->getPoints());}bool GeoSpatialCaculation:: IsLinesIntersect_Ignore_TwoSides(const LineString* lineStringU, const LineString* lineStringV){ if(lineStringU->getPointsCount() <= 1) return false; for (size_t j = 1; j < lineStringU->getPointsCount(); j++) { if(IsLinesIntersect_Ignore_TwoSides(lineStringU->getPoint(j-1), lineStringU->getPoint(j), lineStringV)) return true; } return false;}bool GeoSpatialCaculation:: IsLinesIntersect_Seg_Line(const Point& u1, const Point& u2, const Point& v1, const Point& v2){ return multx(u1,v2,v1) * multx(v2,u2,v1) >= 0;}bool GeoSpatialCaculation:: IsLinesIntersect_Line_Line(const Point& u1, const Point& u2, const Point& v1, const Point& v2){ if(IsPointsInOneLine(u1, u2, v1) && IsPointsInOneLine(u1, u2, v2)) return true; return !IsLinesParallel(u1, u2, v1, v2);}int GeoSpatialCaculation:: PointInLine(const Point& point, const Point& point1, const Point& point2){ double r = multx(point2, point, point1); if(r > ConstValue::ZERO ) return -1; if(r < -ConstValue::ZERO ) return 1; if((point1.X() - point.X()) * (point2.X() - point.X()) <= 0 && (point1.Y() - point.Y()) * (point2.Y() - point.Y()) <= 0) return 0; if(ptDist(point, point2, false) < ptDist(point, point1, false)) return 2; return 3;}int GeoSpatialCaculation:: PointInRectangle(const Point& point, const Envelope& envelope){ if(!envelope.contain(point)) return -1; return point.X() > envelope.getMinX() && point.X() < envelope.getMaxX() && point.Y() > envelope.getMinY() && point.Y() < envelope.getMaxY();}int GeoSpatialCaculation:: PointInCircle(const Point& point, const Point& pointCenter, double radias){ double d = ptDist(pointCenter, point, false); radias *= radias; if(d == radias) return 0; else if(d < radias) return 1; else return -1;}int GeoSpatialCaculation:: PointInEllipse(const Point& point, const Point& pointCenter, double radiasA, double radiasB){ double dx = (point.getX() - pointCenter.getX()) * (point.getX() - pointCenter.getX()); double dy = (point.getY() - pointCenter.getY()) * (point.getY() - pointCenter.getY()); double d = dx / (radiasA*radiasA) + dy / (radiasB*radiasB); if(d == 1) return 0; else if(d < 1) return 1; else return -1;}int GeoSpatialCaculation:: PointInLinearRing(const Point& point, const LinearRing* linearRing){ if(!linearRing->getEnvelope().contain(point)) return -1; if(linearRing->containNode(point)) return 0; if(linearRing->containPoint(point)) return 1; for(int i = 1; i < linearRing->getPointsCount(); i++) { if(IsPointOnLine(point, linearRing->getPoint(i-1), linearRing->getPoint(i))) return 0; } return -1;}int GeoSpatialCaculation:: PointInPolygon(const Point& point, const Polygon* polygon){ if(!polygon->getEnvelope().contain(point)) return -1; if(polygon->containNode(point)) return 0; if(polygon->containPoint(point)) return 1; int n = (int)polygon->getInteriorRingsCount(); for (int j = -1; j < n; j++) { const LinearRing* line = j != -1 ? polygon->getInteriorRing(j) : polygon->getExteriorRing(); for(int i = 1; i < line->getPointsCount(); i++) { if(IsPointOnLine(point, line->getPoint(i-1), line->getPoint(i))) return 0; } } return -1;}bool GeoSpatialCaculation:: IsLineIntersectRectangle(const Point& point1, const Point& point2, const Envelope& envelope){ if(!MakeEnvelope(point1, point2).intersect(envelope)) return false; if(PointInRectangle(point1, envelope) >= 0) return true; if(PointInRectangle(point2, envelope) >= 0) return true; bool bIntersect = IsLinesIntersect(point1, point2, Point(envelope.getMinX(),envelope.getMaxY()), Point(envelope.getMaxX(),envelope.getMinY())); if(bIntersect) return true; bIntersect = IsLinesIntersect(point1, point2, Point(envelope.getMaxX(),envelope.getMaxY()), Point(envelope.getMinX(),envelope.getMinY())); if(bIntersect) return true; return false;}bool GeoSpatialCaculation:: IsLineIntersectCircle(const Point& point1, const Point& point2, const Point& pointCenter, double radias){ if(PointInCircle(point1, pointCenter, radias) >= 0) return true; if(PointInCircle(point2, pointCenter, radias) >= 0) return true; return PtDistToSegment(pointCenter, point1, point2) <= radias;}bool GeoSpatialCaculation:: IsLineIntersectEllipse(const Point& point1, const Point& point2, const Point& pointCenter, double radiasA, double radiasB){ if(!MakeEnvelope(point1, point2).intersect(MakeEnvelope(pointCenter, radiasA > radiasB ? radiasA : radiasB))) return false; if(PointInEllipse(point1, pointCenter, radiasA, radiasB) >= 0) return true; if(PointInEllipse(point2, pointCenter, radiasA, radiasB) >= 0) return true; PointVec pointVecIntersect; lineIntersectEllipseP(point1, point2, pointCenter, radiasA, radiasB, pointVecIntersect); return pointVecIntersect.size() > 0;}bool GeoSpatialCaculation:: IsLineIntersectPolygon(const Point& point1, const Point& point2, const Polygon* polygon){ if(!MakeEnvelope(point1, point2).intersect(polygon->getEnvelope())) return false; if(PointInPolygon(point1, polygon) >= 0) return true; if(PointInPolygon(point2, polygon) >= 0) return true; if(IsLinesIntersect(point1, point2, polygon->getExteriorRing())) return true; for (size_t i = 0; i < polygon->getInteriorRingsCount(); i++) { if(IsLinesIntersect(point1, point2, polygon->getInteriorRing(i))) return true; } return false; }bool GeoSpatialCaculation:: IsLineIntersectRectangle(const LineString* lineString, const Envelope& envelope){ int n = (int)lineString->getPointsCount(); if(n <= 1) return false; for(int i = 1; i < n; i++) { if(IsLineIntersectRectangle(lineString->getPoint(i-1), lineString->getPoint(i), envelope)) return true; } return false;}bool GeoSpatialCaculation:: IsLineIntersectCircle(const LineString* lineString, const Point& pointCenter, double radias){ int n = (int)lineString->getPointsCount(); if(n <= 1) return false; for(int i = 1; i < n; i++) { if(IsLineIntersectCircle(lineString->getPoint(i-1), lineString->getPoint(i), pointCenter, radias)) return true; } return false;}bool GeoSpatialCaculation:: IsLineIntersectEllipse(const LineString* lineString, const Point& pointCenter, double radiasA, double radiasB){ int n = (int)lineString->getPointsCount(); if(n <= 1) return false; for(int i = 1; i < n; i++) { if(IsLineIntersectEllipse(lineString->getPoint(i-1), lineString->getPoint(i), pointCenter, radiasA, radiasB)) return true; } return false;}bool GeoSpatialCaculation:: IsLineIntersectPolygon(const LineString* lineString, const Polygon* polygon){ int n = (int)lineString->getPointsCount(); if(n <= 1) return false; for(int i = 1; i < n; i++) { if(IsLineIntersectPolygon(lineString->getPoint(i-1), lineString->getPoint(i), polygon)) return true; } return false;}bool GeoSpatialCaculation:: IsLineInRectangle(const Point& point1, const Point& point2, const Envelope& envelope){ if(PointInRectangle(point1, envelope) == -1) return false; return PointInRectangle(point2, envelope) != -1;}bool GeoSpatialCaculation:: IsLineInCircle(const Point& point1, const Point& point2, const Point& pointCenter, double radias){ if(PointInCircle(point1, pointCenter, radias) == -1) return false; return PointInCircle(point2, pointCenter, radias) != -1;}bool GeoSpatialCaculation:: IsLineInEllipse(const Point& point1, const Point& point2, const Point& pointCenter, double radiasA, double radiasB){ if(PointInEllipse(point1, pointCenter, radiasA, radiasB) == -1) return false; return PointInEllipse(point2, pointCenter, radiasA, radiasB) != -1;}bool GeoSpatialCaculation:: IsLineInPolygon(const Point& point1, const Point& point2, const Polygon* polygon){ Envelope ev = MakeEnvelope(point1, point2); if(!polygon->getEnvelope().intersect(ev)) return false; if(PointInPolygon(point1, polygon) == -1 || PointInPolygon(point2, polygon) == -1) return false; int n = (int)polygon->getInteriorRingsCount(); for(int i = -1; i < n; i++) { const LinearRing* line = (i == -1) ? polygon->getExteriorRing() : polygon->getInteriorRing(i); if(line == 0) continue; if(!line->getEnvelope().contain(ev)) continue; for(size_t j = 1; j < line->getPointsCount(); j++) { PointVec pointVecIntersect; bool bIntersect = GetIntersectPoint(point1, point2, line->getPoint(j-1), line->getPoint(j), pointVecIntersect); if(!bIntersect) continue; Point pt; for(size_t k = 0; k < pointVecIntersect.size(); k++) { Point p = pointVecIntersect[k]; Point pt = p == point1 ? point2 : point1; Point p1 = PointExtend(pt, p, ConstValue::ZERO * 2); if(IsPointOnLine(p1, point1, point2) && PointInPolygon(p1, polygon) == -1) return false; Point p2 = PointExtend(pt, p, -ConstValue::ZERO * 2); if(IsPointOnLine(p2, point1, point2) && PointInPolygon(p2, polygon) == -1) return false; } } } return true;}bool GeoSpatialCaculation:: IsLineInRectangle(const LineString* lineString, const Envelope& envelope){ int n = (int)lineString->getPointsCount(); if(n <= 1) return false; for(int i = 1; i < n; i++) { if(IsLineInRectangle(lineString->getPoint(i-1), lineString->getPoint(i), envelope)) return true; } return false;}bool GeoSpatialCaculation:: IsLineInCircle(const LineString* lineString, const Point& pointCenter, double radias){ int n = (int)lineString->getPointsCount(); if(n <= 1) return false; for(int i = 1; i < n; i++) { if(IsLineInCircle(lineString->getPoint(i-1), lineString->getPoint(i), pointCenter, radias)) return true; } return false;}bool GeoSpatialCaculation:: IsLineInEllipse(const LineString* lineString, const Point& pointCenter, double radiasA, double radiasB){ int n = (int)lineString->getPointsCount(); if(n <= 1) return false; for(int i = 1; i < n; i++) { if(IsLineInEllipse(lineString->getPoint(i-1), lineString->getPoint(i), pointCenter, radiasA, radiasB)) return true; } return false;}bool GeoSpatialCaculation:: IsLineInPolygon(const LineString* lineString, const Polygon* polygon){ int n = (int)lineString->getPointsCount(); if(n <= 1) return false; for(int i = 1; i < n; i++) { if(IsLineInPolygon(lineString->getPoint(i-1), lineString->getPoint(i), polygon)) return true; } return false;}bool GeoSpatialCaculation:: IsRectangleIntersectRectangle(const Envelope& envelope1, const Envelope& envelope2){ return envelope1.intersect(envelope2);}bool GeoSpatialCaculation:: IsRectangleIntersectCircle(const Envelope& envelope, const Point& pointCenter, double radias){ if(!envelope.intersect(MakeEnvelope(pointCenter, radias))) return false; if(PointInRectangle(pointCenter, envelope) != -1) return true; return PtDistToRectangle(pointCenter, envelope, false) <= radias * radias;}bool GeoSpatialCaculation:: IsRectangleIntersectPolygon(const Envelope& envelope, const Polygon* polygon){ if(!envelope.intersect(polygon->getEnvelope())) return false; int n = (int)polygon->getInteriorRingsCount(); for(int i = -1; i < n; i++) { const LinearRing* line = (i == -1) ? polygon->getExteriorRing() : polygon->getInteriorRing(i); if(line == 0) continue; if(!envelope.intersect(line->getEnvelope())) continue; for(size_t j = 1; j < line->getPointsCount(); j++) { if(IsLineIntersectRectangle(line->getPoint(j-1), line->getPoint(j), envelope)) return true; } } return false;}bool GeoSpatialCaculation:: IsCircleIntersectCircle(const Point& pointCenter1, double radias1, const Point& pointCenter2, double radias2){ return ptDist(pointCenter1, pointCenter2) <= radias1 + radias2;}bool GeoSpatialCaculation:: IsCircleIntersectPolygon(const Point& pointCenter, double radias, const Polygon* polygon){ if(!polygon->getEnvelope().intersect(MakeEnvelope(pointCenter, radias))) return false; if(PointInPolygon(pointCenter, polygon) >= 0) return true; return PtDistToPolygon(pointCenter, polygon) <= radias;}bool GeoSpatialCaculation:: IsPolygonIntersectPolygon(const Polygon* polygon1, const Polygon* polygon2){ Envelope evPoly1 = polygon1->getEnvelope(); Envelope evPoly2 = polygon2->getEnvelope(); if(!evPoly1.intersect(evPoly2)) return false; const Polygon* smallPoly = evPoly2.getAera() > evPoly1.getAera() ? polygon1 : polygon2; const Polygon* largePoly = smallPoly == polygon1 ? polygon2 : polygon1; int n = (int)smallPoly->getInteriorRingsCount(); for(int i = -1; i < n; i++) { const LinearRing* line1 = (i == -1) ? smallPoly->getExteriorRing() : smallPoly->getInteriorRing(i); if(line1 == 0) continue; Envelope evLine1 = line1->getEnvelope(); if(!IsRectangleIntersectRectangle(evLine1, evPoly2)) continue; for(size_t j = 1; j < line1->getPointsCount(); j++) { if(IsLineIntersectPolygon(line1->getPoint(j-1), line1->getPoint(j), largePoly)) return true; } } return false;}bool GeoSpatialCaculation:: IsRectangleInRectangle(const Envelope& envelope1, const Envelope& envelope2){ return envelope2.contain(envelope1);}bool GeoSpatialCaculation:: IsRectangleInCircle(const Envelope& envelope, const Point& pointCenter, double radias){ PointVec points = getEnvelopePoints(envelope); for (size_t i = 0; i < 4; i++) { if(PointInCircle(points[i], pointCenter, radias) == -1) return false; } return true;}bool GeoSpatialCaculation:: IsRectangleInPolygon(const Envelope& envelope, const Polygon* polygon){ if(!polygon->getEnvelope().contain(envelope)) return false; if(!IsLineInPolygon(Point(envelope.getMinX(), envelope.getMaxY()), Point(envelope.getMaxX(), envelope.getMaxY()), polygon)) return false; if(!IsLineInPolygon(Point(envelope.getMaxX(), envelope.getMaxY()), Point(envelope.getMaxX(), envelope.getMinY()), polygon)) return false; if(!IsLineInPolygon(Point(envelope.getMaxX(), envelope.getMinY()), Point(envelope.getMinX(), envelope.getMinY()), polygon)) return false; if(!IsLineInPolygon(Point(envelope.getMinX(), envelope.getMinY()), Point(envelope.getMinX(), envelope.getMaxY()), polygon)) return false; return true;}bool GeoSpatialCaculation:: IsCircleInRectangle(const Point& pointCenter, double radias, const Envelope& envelope){ if(PointInRectangle(pointCenter, envelope) != 1) return false; return -PtDistToRectangle(pointCenter, envelope) >= radias;}bool GeoSpatialCaculation:: IsCircleInCircle(const Point& pointCenter1, double radias1, const Point& pointCenter2, double radias2){ return ptDist(pointCenter1, pointCenter2) + radias1 <= radias2;}bool GeoSpatialCaculation:: IsCircleInPolygon(const Point& pointCenter, double radias, const Polygon* polygon){ if(PointInPolygon(pointCenter, polygon) != 1) return false; return -PtDistToPolygon(pointCenter, polygon) >= radias;}bool GeoSpatialCaculation:: IsPolygonInRectangle(const Polygon* polygon, const Envelope& envelope){ if(!envelope.contain(polygon->getEnvelope())) return false; const LinearRing* line = polygon->getExteriorRing(); if(line == 0 || line->getPointsCount() == 0) return false; for(size_t j = 0; j < line->getPointsCount(); j++) { if(PointInRectangle(line->getPoint(j), envelope) == -1) return false; } return true;}bool GeoSpatialCaculation:: IsPolygonInCircle(const Polygon* polygon, const Point& pointCenter, double radias){ if(!MakeEnvelope(pointCenter, radias).contain(polygon->getEnvelope())) return false; const LinearRing* line = polygon->getExteriorRing(); if(line == 0 || line->getPointsCount() == 0) return false; for(size_t j = 0; j < line->getPointsCount(); j++) { if(PointInCircle(line->getPoint(j), pointCenter, radias) == -1) return false; } return true;}bool GeoSpatialCaculation:: IsPolygonInPolygon(const Polygon* polygon1, const Polygon* polygon2){ if(!polygon2->getEnvelope().contain(polygon1->getEnvelope())) return false; const LinearRing* line = polygon1->getExteriorRing(); if(line == 0 || line->getPointsCount() == 0) return false; for(size_t j = 1; j < line->getPointsCount(); j++) { if(!IsLineInPolygon(line->getPoint(j-1), line->getPoint(j), polygon2)) return false; } return true;}//***************************** 3.求值 ********************************//double GeoSpatialCaculation:: PointPrj2Line(double x,double y,double x1,double y1,double x2,double y2,double &XResult,double &YResult){ double dis = ptDist(x1,y1, x2,y2, false); double r = 0; if(dis != 0) r = multd(Point(x,y), Point(x2,y2), Point(x1,y1)) / dis; XResult = x1 + r*(x2 - x1); YResult = y1 + r*(y2 - y1); return r;}double GeoSpatialCaculation:: PointPrj2Line(const Point& point, const Point& point1, const Point& point2, Point& pointRet){ double XResult = 0, YResult = 0; double r = PointPrj2Line( point.getX(), point.getY(), point1.getX(), point1.getY(), point2.getX(), point2.getY(), XResult, YResult ); pointRet.setX(XResult); pointRet.setY(YResult); return r;}void GeoSpatialCaculation:: PointRotate(double xLoc, double yLoc, double xOld, double yOld, double& xNew, double& yNew, double angle){ xNew = xLoc + cos(angle) * (xOld - xLoc) - sin(angle) * (yOld - yOld); yNew = yLoc + sin(angle) * (xOld - xLoc) + cos(angle) * (yOld - yOld);}Point GeoSpatialCaculation:: PointRotate(const Point& pointLoc, const Point& pointOld, double angle){ double x = pointLoc.X() + cos(angle) * (pointOld.X() - pointLoc.X()) - sin(angle) * (pointOld.Y() - pointLoc.Y()); double y = pointLoc.Y() + sin(angle) * (pointOld.X() - pointLoc.X()) + cos(angle) * (pointOld.Y() - pointLoc.Y()); return Point(x, y);}void GeoSpatialCaculation:: PointOnLine(double x1, double y1, double x2, double y2, double dScale,double &XResult,double &YResult){ if (fabs(x1 - x2) < ConstValue::ZERO) { XResult = x1 ; YResult = y1 + (y2 - y1) * dScale; } else { XResult = x1 + (x2 - x1) * dScale; YResult = y1 + (y2 - y1) * dScale; }}Point GeoSpatialCaculation:: PointOnLine(const Point& point1, const Point& point2, double dScale){ double XResult = 0, YResult = 0; PointOnLine( point1.getX(), point1.getY(), point2.getX(), point2.getY(), dScale, XResult, YResult); return Point(XResult, YResult);}Point GeoSpatialCaculation:: PointOnLine(const LineString* lineString, double dScale){ return PointOnLine(lineString->getPoints(), dScale);}Point GeoSpatialCaculation:: PointOnLine(const PointVec& pointVec, double dScale){ Point resultPoint; if(dScale == 0) return Point(pointVec[0]); else if(dScale == 1) return Point(pointVec[pointVec.size()-1]); double totalLength = LineLength(pointVec); double needLength = totalLength * dScale; double length = 0; size_t pointsCount = pointVec.size() - 1; for (size_t i = 0; i < pointsCount; i++) { const Point& point1 = pointVec[i]; const Point& point2 = pointVec[i+1]; length += ptDist(point1, point2); if(length >= needLength) { if(length == needLength) resultPoint = Point(point2); else { length -= needLength; length /= ptDist(point1, point2); resultPoint = PointOnLine(point1, point2, 1-length); } break; } } return resultPoint;}double GeoSpatialCaculation:: ScaleOnLine(const Point& point, const PointVec& pointVec){ int segIndex = -1; if(!IsPointOnLine(point, pointVec, segIndex)) return -1; if(segIndex < 0 || segIndex > pointVec.size() - 1) return -1; double d = LineLength(pointVec, 0, segIndex+1) + ptDist(point, pointVec[segIndex+1]); return d / LineLength(pointVec);}double GeoSpatialCaculation:: ScaleOnLine(const Point& point, const LineString* lineString){ return ScaleOnLine(point, lineString->getPoints());}void GeoSpatialCaculation:: PointCenter(double x1, double y1, double x2, double y2, double& x, double& y){ ptCenter(x1, y1, x2, y2, x, y);}Point GeoSpatialCaculation:: PointCenter(const Point& point1, const Point& point2){ double x = 0;double y = 0; ptCenter(point1.getX(), point1.getY(), point2.getX(), point2.getY(), x, y); return Point(x, y);}void GeoSpatialCaculation:: PointExtend(double x1, double y1, double x2, double y2, double extendDistance, double &XResult,double &YResult){ double dis = ptDist(x1, y1, x2, y2); if(dis < ConstValue::ZERO) return; double scale = 1 + extendDistance / dis; PointOnLine(x1, y1, x2, y2, scale, XResult, YResult);}Point GeoSpatialCaculation:: PointExtend(const Point& point1, const Point& point2, double extendDistance){ double XResult = 0, YResult = 0; PointExtend( point1.getX(), point1.getY(), point2.getX(), point2.getY(), extendDistance, XResult, YResult); return Point(XResult, YResult);}Point GeoSpatialCaculation:: PointExtend(const LineString* lineString, double extendDistance, bool beginStart){ Point pt(ConstValue::MAXVALUE, ConstValue::MAXVALUE); if(lineString == 0) return pt; return PointExtend(lineString->getPoints(), extendDistance, beginStart);}Point GeoSpatialCaculation:: PointExtend(const PointVec& pointVec, double extendDistance, bool beginStart){ if(extendDistance >= 0) { if(!beginStart) return PointExtend(pointVec[pointVec.size()-2], pointVec[pointVec.size()-1], extendDistance); else return PointExtend(pointVec[1], pointVec[0], extendDistance); } else { double len = LineLength(pointVec); double sacle = -extendDistance / len; if(!beginStart) sacle = 1 - sacle; if(sacle > 1) { return PointExtend(pointVec[pointVec.size()-2], pointVec[pointVec.size()-1], -extendDistance-len); } else if(sacle < 0) { return PointExtend(pointVec[1], pointVec[0], -extendDistance-len); } return PointOnLine(pointVec, sacle); }}void GeoSpatialCaculation:: PPLDN(double x,double y,double x1,double y1,double x2, double y2,double distance, double& resultx,double& resulty){ double a1,a2; a1 = LineOrientation(x1, y1, x2, y2); a2 = a1 + ConstValue::PI / 2 ; if (a2 > 2 * ConstValue::PI) a2 -= 2 * ConstValue::PI ; resultx = x + (fabs(distance) * cos(a2)) ; resulty = y + (fabs(distance) * sin(a2)) ; if (distance < -ConstValue::ZERO) { resultx = 2 * x - resultx ; resulty = 2 * y - resulty ; }}void GeoSpatialCaculation:: PPLDN(const Point& point, const Point& point1, const Point& point2, Point& pointResult, double d){ double XResult = 0, YResult = 0; PPLDN( point.getX(), point.getY(), point1.getX(), point1.getY(), point2.getX(), point2.getY(), d, XResult, YResult ); pointResult.setX(XResult); pointResult.setY(YResult);}void GeoSpatialCaculation:: P2LSymmetry(const Point& point1, const Point& point2, const Point& point3, Point& pointResult){ double xRef = 0, yRef = 0; PointPrj2Line(point3.getX(), point3.getY(), point1.getX(), point1.getY(), point2.getX(), point2.getY(), xRef, yRef); pointResult.setX(2 * xRef - point3.getX()); pointResult.setY(2 * yRef - point3.getY());}void GeoSpatialCaculation:: P2LSymmetry(double x1, double y1, double x2, double y2, double x3, double y3, double &XResult, double &YResult){ double X = 0, Y = 0;//投影点坐标 PointPrj2Line(x3, y3, x1, y1, x2, y2, X, Y); P2PSymmetry(x3, y3, X, Y, XResult, YResult);}void GeoSpatialCaculation:: P2PSymmetry(double x,double y,double symmetryX,double symmetryY, double& resultx,double& resulty){ resultx = 2 * symmetryX - x; resulty = 2 * symmetryY - y;}void GeoSpatialCaculation:: P2PSymmetry(const Point& point, const Point& symmetryPoint, Point& pointResult){ pointResult.setX(2 * symmetryPoint.getX() - point.getX()); pointResult.setY(2 * symmetryPoint.getY() - point.getY());}bool GeoSpatialCaculation:: GetIntersectPoint(const Point& u1, const Point& u2, const Point& v1, const Point& v2, PointVec& pointVecIntersect){ pointVecIntersect.clear(); if(!IsLinesIntersect(u1, u2, v1, v2)) return false; if(IsLinesParallel(u1, u2, v1, v2)) //如果两线平行且重合，取重合部分的两个端点 { if(IsPointOnLine(u1, v1, v2)) pointVecIntersect.push_back(u1); if(IsPointOnLine(u2, v1, v2)) pointVecIntersect.push_back(u2); if(pointVecIntersect.size() == 2) return true; if(v1 != u1 && v1 != u2 && IsPointOnLine(v1, u1, u2)) pointVecIntersect.push_back(v1); if(pointVecIntersect.size() == 2) return true; if(v2 != u1 && v2 != u2 && IsPointOnLine(v2, u1, u2)) pointVecIntersect.push_back(v2); return pointVecIntersect.size() > 0; } Point pt; bool bIntersect = getIntersectPoint(u1, u2, v1, v2, pt); if(bIntersect) pointVecIntersect.push_back(pt); return bIntersect;}bool GeoSpatialCaculation:: GetIntersectPoint(const Point& u1, const Point& u2, const LineString* lineString, PointVec& pointVecIntersect){ pointVecIntersect.clear(); if(lineString == 0 || lineString->getPointsCount() <= 1) return false; for (size_t j = 1; j < lineString->getPointsCount(); j++) { PointVec ptsRet; if(GetIntersectPoint(u1, u2, lineString->getPoint(j-1), lineString->getPoint(j), ptsRet)) pointVecIntersect.insert(pointVecIntersect.end(), ptsRet.begin(), ptsRet.end()); } uniquePoints(pointVecIntersect); return pointVecIntersect.size() != 0;}bool GeoSpatialCaculation:: GetIntersectPoint_Seg_Line(const Point& u1, const Point& u2, const Point& v1, const Point& v2, PointVec& pointVecIntersect){ pointVecIntersect.clear(); if(!IsLinesIntersect_Seg_Line(u1, u2, v1, v2)) return false; if(IsLinesParallel(u1, u2, v1, v2)) { pointVecIntersect.push_back(u1); pointVecIntersect.push_back(u2); return true; } Point pt; bool bIntersect = getIntersectPoint(u1, u2, v1, v2, pt); if(bIntersect) pointVecIntersect.push_back(pt); return bIntersect;}bool GeoSpatialCaculation:: GetIntersectPoint_Line_Line(const Point& u1, const Point& u2, const Point& v1, const Point& v2, Point& pointRet){ if(IsLinesParallel(u1, u2, v1, v2)) return false; return getIntersectPoint(u1, u2, v1, v2, pointRet);}void GeoSpatialCaculation:: GetIntersectPoint_Seg_Rectangle(const Point& point1, const Point& point2, const Envelope& envelope, PointVec& pointVecIntersect){ pointVecIntersect.clear(); if(!MakeEnvelope(point1, point2).intersect(envelope)) return; int n1 = PointInRectangle(point1, envelope); int n2 = PointInRectangle(point2, envelope); if(n1 == 1 && n2 == 1) return; PointVec points = getEnvelopePoints(envelope); for (size_t i = 0; i < 4; i++) { PointVec pts; bool bIntersect = false; if(i != 3) bIntersect = GetIntersectPoint(point1, point2, points[i], points[i+1], pts); else bIntersect = GetIntersectPoint(point1, point2, points[3], points[0], pts); if(bIntersect) pointVecIntersect.insert(pointVecIntersect.end(), pts.begin(), pts.end()); } uniquePoints(pointVecIntersect);}void GeoSpatialCaculation:: GetIntersectPoint_Seg_Circle(const Point& point1, const Point& point2, const Point& pointCenter, double radias, PointVec& pointVecIntersect){ pointVecIntersect.clear(); Point pt; PointPrj2Line(pointCenter, point1, point2, pt); double dis2 = ptDist(pointCenter, pt, false); double radias2 = radias * radias; if(dis2 > radias2) return; double d = -sqrt(radias2 - dis2); Point ps = point1 == pt ? PointCenter(point1, point2) : point1; Point pt1 = PointExtend(ps, pt, d); ps = point2 == pt ? PointCenter(point1, point2) : point2; Point pt2 = PointExtend(ps, pt, d); if(IsPointOnLine(pt1, point1, point2)) pointVecIntersect.push_back(pt1); if(IsPointOnLine(pt2, point1, point2) && pt1 != pt2) pointVecIntersect.push_back(pt2);}void GeoSpatialCaculation:: GetIntersectPoint_Seg_Ellipse(const Point& point1, const Point& point2, const Point& pointCenter, double radiasA, double radiasB, PointVec& pointVecIntersect){ pointVecIntersect.clear(); int n1 = PointInEllipse(point1, pointCenter, radiasA, radiasB); int n2 = PointInEllipse(point2, pointCenter, radiasA, radiasB); if(n1 == 1 && n2 == 1) return; if(n1 == 0 && n2 == 0) { pointVecIntersect.push_back(point1); pointVecIntersect.push_back(point2); } else if(n1 == 0 || n2 == 0) pointVecIntersect.push_back(n1 == 0 ? point1 : point2); else if((n1 == 1 && n2 == -1) || (n1 == -1 && n2 == 1)) lineIntersectEllipseP(point1, point2, pointCenter, radiasA, radiasB, pointVecIntersect); else { Point pt; double dis = PtDistToSegment(pointCenter, point1, point2, pt); int n = PointInEllipse(pt, pointCenter, radiasA, radiasB); if(n == 0) pointVecIntersect.push_back(pt); else if(n == 1) lineIntersectEllipseP(point1, point2, pointCenter, radiasA, radiasB, pointVecIntersect); }}void GeoSpatialCaculation:: GetIntersectPoint_Seg_Polygon(const Point& point1, const Point& point2, const Polygon* polygon, PointVec& pointVecIntersect){ PointVec linearingIndexAndSegIndex; getIntersectPoint_Seg_Polygon(point1, point2, polygon, pointVecIntersect, linearingIndexAndSegIndex, false);}void GeoSpatialCaculation:: GetIntersectPoint_Seg_Polygon(const Point& point1, const Point& point2, const Polygon* polygon, PointVec& pointVecIntersect, PointVec& linearingIndexAndSegIndex){ getIntersectPoint_Seg_Polygon(point1, point2, polygon, pointVecIntersect, linearingIndexAndSegIndex, true);}int GeoSpatialCaculation:: PointConvexity(Point pointCur, Point ptFirst, Point ptLast, bool isClockWise){ if(IsPointsEqual(ptLast, pointCur) || IsPointsEqual(ptLast, ptFirst) || IsPointsEqual(pointCur, ptFirst)) return 0; double dmultx = multx(pointCur, ptLast, ptFirst) ; if(fabs(dmultx) <= ConstValue::ZERO) return 0; if((isClockWise && dmultx > 0) || (!isClockWise && dmultx < 0)) return -1;//凹部 return 1;}int GeoSpatialCaculation:: PointConvexity(int pointIndex, const LinearRing* linearRing){ if(pointIndex == 0 || pointIndex == linearRing->getPointsCount()-1) return PointConvexity(linearRing->getPoint(0), linearRing->getPoint(linearRing->getPointsCount()-2), linearRing->getPoint(1), linearRing->isClockwise()); return PointConvexity(linearRing->getPoint(pointIndex), linearRing->getPoint(pointIndex-1), linearRing->getPoint(pointIndex+1), linearRing->isClockwise());}int GeoSpatialCaculation:: PointConvexity(int pointIndex, const Polygon* polygon, int lineRingIndex){ return PointConvexity(pointIndex, lineRingIndex == -1 ? polygon->getExteriorRing() : polygon->getInteriorRing(lineRingIndex));}int GeoSpatialCaculation:: PointConvexity(Point pointCur, Point ptFirst, Point ptLast){ if(IsPointsEqual(ptLast, pointCur) || IsPointsEqual(ptLast, ptFirst) || IsPointsEqual(pointCur, ptFirst)) return 0; double dmultx = multx(pointCur, ptLast, ptFirst) ; if(fabs(dmultx) <= ConstValue::ZERO) return 0; if(dmultx < 0) return 1; else return -1;}int GeoSpatialCaculation:: PointConvexity(int pointIndex, const LineString* lineString){ if(lineString->getGeometryType() == GeometryType_LinearRing) return PointConvexity(pointIndex, (const LinearRing*)lineString); if(pointIndex <= 0 || pointIndex >= lineString->getPointsCount() -1) return 0; return PointConvexity(lineString->getPoint(pointIndex), lineString->getPoint(pointIndex-1), lineString->getPoint(pointIndex+1));}void GeoSpatialCaculation:: PointsConvexHull(const PointVec& points, PointVec& pointVecRet){ PointVec pointsRet; pointsRet.insert(pointsRet.begin(), points.begin(), points.end()); uniquePoints(pointsRet); if(pointsRet.size() <= 3) { pointVecRet.assign(points.begin(), points.end()); return; } std::vector stack; stack.push_back(0); stack.push_back(1); int top=2; int n = pointsRet.size(); for (int i=2;i 1 && multx(pointsRet[i],pointsRet[stack[top-1]],pointsRet[stack[top-2]])>0) top--; if(stack.size() <= top) stack.push_back(i); else stack[top] = i; top++; } int t=top; //重新赋值 if(stack.size() <= top) stack.push_back(n-2); else stack[top] = n-2; top++; for (int i = n-3;i>=0;i--)//搜索最后三个点 { while (top>t && multx(pointsRet[i],pointsRet[stack[top-1]],pointsRet[stack[top-2]])>0) top--; if(stack.size() <= top) stack.push_back(i); else stack[top] = i; top++; } for (int i=0;i getExteriorRing(); if(linearRing == 0 || linearRing->getPointsCount() < 3) return pt; bool isClockwise = linearRing->isClockwise(); for(size_t i = 2; i < linearRing->getPointsCount(); i++) { Point pt0 = linearRing->getPoint(i-2); Point pt1 = linearRing->getPoint(i-1); Point pt2 = linearRing->getPoint(i); if(PointConvexity(pt1, pt0, pt2, isClockwise) <= 0) continue;//共线或者凹部 Point center = PointCenter(pt0, pt2); Point ptCur = PointCenter(pt1, center); PointVec ptVec; GetIntersectPoint_Seg_Polygon(pt1, ptCur, polygon, ptVec); if(ptVec.size() <= 2) return center; Point ptRet = getNearistPoint(pt1, ptVec, true); if(ptRet.getX() != ConstValue::MAXVALUE) return PointCenter(pt1, ptRet); } return pt;}Point GeoSpatialCaculation:: LinearRingCentroid(const LinearRing* linearRing){ PointVec p = linearRing->getPoints(); double area = 0, cx = 0, cy = 0; for(size_t i = 0; i < p.size()-1; i++) { double d = p[i].getX()*p[i+1].getY() - p[i+1].getX()*p[i].getY(); area += d * 0.5; cx += d * (p[i].getX() + p[i+1].getX()); cy += d * (p[i].getY() + p[i+1].getY()); } return Point(cx /= 6*area, cy /= 6*area);}Point GeoSpatialCaculation:: PolygonCentroid(const Polygon* polygon){ Point pt = LinearRingCentroid(polygon->getExteriorRing()); if(PointInPolygon(pt, polygon) != 1) { double radias = 0; return PolygonMaxInnerCircle(polygon, radias); } return pt;}Point GeoSpatialCaculation:: LinearRingMaxInnerCircle(const LinearRing* linearRing, double& radias, int iterateCount){ Polygon* polygon = GeometryFactory::createPolygon(linearRing, false); Point point = PolygonMaxInnerCircle(polygon, radias, iterateCount); GeometryFactory::destoryGeometry(polygon); return point;}Point GeoSpatialCaculation:: PolygonMaxInnerCircle(const Polygon* polygon, double& radias, int iterateCount){ radias = 0; Point ptRet(ConstValue::MAXVALUE, ConstValue::MAXVALUE); const LinearRing* line = polygon->getExteriorRing(); if(line == 0 || line->getPointsCount() < 3) return ptRet; //随机点求最大内圆 PointVec points = getCentralPoints(line->getEnvelope().getMinX(), line->getEnvelope().getMinY(), line->getEnvelope().getMaxX(), line->getEnvelope().getMaxY(), iterateCount); double radias2 = ConstValue::MAXVALUE; bool bFirstPoint = true; for(size_t i = 0; i < points.size(); i++) { Point point = points[i]; if(PointInPolygon(point, polygon) != 1) continue; double dmin = ConstValue::MAXVALUE; Point ptMin ; bool bsus = true; for (int j = -1; j < (int)polygon->getInteriorRingsCount(); j++) { const LinearRing* linearRing = (j == -1) ? polygon->getExteriorRing() : polygon->getInteriorRing(j); for (size_t k = 1; k < linearRing->getPointsCount(); k++) { double dis = PtDistToSegment(point, linearRing->getPoint(k-1), linearRing->getPoint(k), false); if(!bFirstPoint && dis < radias2) { bsus = false; break; } if(bFirstPoint) { if(radias2 > dis) { radias2 = dis; ptRet = point; } } else { if(dmin > dis) { dmin = dis; ptMin = point; } } } } if(bsus && !bFirstPoint && radias2 < dmin) { radias2 = dmin; ptRet = ptMin; } if(bFirstPoint) bFirstPoint = false; } if(radias2 == ConstValue::MAXVALUE) { ptRet = PolygonInnerPoint(polygon); radias2 = PtDistToPolygon(ptRet, polygon, false); } radias = sqrt(radias2); return ptRet;}double GeoSpatialCaculation:: LineOrientation(double x1, double y1, double x2, double y2){ double dx = x2 - x1; double dy = y2 - y1; if (dx < ConstValue::ZERO && dy < ConstValue::ZERO) { return 0; } double magnitude = dx * dx + dy * dy; double cosA = dx / magnitude; double arcCosA = Mathematics::Math::ArcCos(cosA); return y2 < y1 ? ConstValue::PI * 2 - arcCosA : arcCosA;}double GeoSpatialCaculation:: LineOrientation(const Point& point1, const Point& point2){ return LineOrientation(point1.getX(), point1.getY(), point2.getX(), point2.getY());}double GeoSpatialCaculation:: LinesAngle(const Point& anglePoint, const Point& point1, const Point& point2){ if(point1 == anglePoint || point2 == anglePoint) return 0; return Nmo::Geometries::Vector::Angle(point1 - anglePoint, point2 - anglePoint);}double GeoSpatialCaculation:: LinesAngle(const Point& startPoint1, const Point& endPoint1, const Point& startPoint2, const Point& endPoint2){ if(startPoint1 == endPoint1 || startPoint2 == endPoint2) return 0; return Nmo::Geometries::Vector::Angle(startPoint1 - endPoint1, startPoint2 - endPoint2);}Point GeoSpatialCaculation:: LinesHalfAnglePoint(const Point& anglePoint, const Point& point1, const Point& point2, bool clockWise){ double or1 = LineOrientation(anglePoint, point1); double or2 = LineOrientation(anglePoint, point2); double angle = (or2 - or1) * 0.5; if(clockWise) { if(angle > 0) angle -= ConstValue::PI; } else { if(angle < 0) angle += ConstValue::PI; } return PointRotate(anglePoint, point1, angle);}double GeoSpatialCaculation:: LineLength(const LineString* lineString, int ptFromIndex, int ptToIndex){ if(lineString == 0) return 0; PointVec pointVec = lineString->getPoints(); return LineLength(pointVec, ptFromIndex, ptToIndex);}double GeoSpatialCaculation:: LineLength(const PointVec& pointVec, int ptFromIndex, int ptToIndex){ if(pointVec.size() <= 1) return 0; if(ptFromIndex < 0) ptFromIndex = 0; if(ptToIndex >= pointVec.size()) ptToIndex = pointVec.size() - 1; if(ptFromIndex >= ptToIndex) return 0; double len = 0; for(int i = ptFromIndex; i < ptToIndex; i++) len += ptDist(pointVec[i], pointVec[i+1]); return len;}double GeoSpatialCaculation:: PolygonLength(const Polygon* polygon){ if(polygon == 0 || polygon->getExteriorRing() == 0) return 0; const LinearRing* linearRing = polygon->getExteriorRing(); double len = linearRing->getLength(); for (size_t j = 0; j < polygon->getInteriorRingsCount(); j++) len += polygon->getInteriorRing(j) ? polygon->getInteriorRing(j)->getLength() : 0; return len;}double GeoSpatialCaculation:: MultiPolygonLength(const MultiPolygon* multiPolygon){ if(multiPolygon == 0) return 0; double len = 0; for (size_t j = 0; j < multiPolygon->getCount(); j++) len += PolygonLength(multiPolygon->getPolygon(j)); return len;}double GeoSpatialCaculation:: MultiPolygonArea(const MultiPolygon* multiPolygon){ double area = 0; for (size_t i = 0; i < multiPolygon->getCount(); ++i) area += multiPolygon->getPolygon(i)->getAera(); return area;}double GeoSpatialCaculation:: TriangleAera(const Point& point1, const Point& point2, const Point& point3){ double AeraTimes = Vector::Cross(Vector(point1.getX(), point1.getY()), Vector(point2.getX(), point2.getY())); AeraTimes += Vector::Cross(Vector(point2.getX(), point2.getY()), Vector(point3.getX(), point3.getY())); AeraTimes += Vector::Cross(Vector(point3.getX(), point3.getY()), Vector(point1.getX(), point1.getY())); AeraTimes = AeraTimes >= 0 ? AeraTimes : -AeraTimes; return AeraTimes * 0.5;}double GeoSpatialCaculation:: PointVecAera(const PointVec& pointVec){ if(pointVec.size() < 3) return 0.0; double AeraTimes = 0.0; for (size_t i = 0; i < pointVec.size(); ++i) { const Point& FrontPoint = pointVec[i]; const Point& BackPoint = pointVec[(i + 1) % pointVec.size()]; AeraTimes += Vector::Cross(Vector(FrontPoint.getX(), FrontPoint.getY()), Vector(BackPoint.getX(), BackPoint.getY())); } AeraTimes = AeraTimes >= 0 ? AeraTimes : -AeraTimes; return AeraTimes * 0.5;}Envelope GeoSpatialCaculation:: MakeEnvelope(const Point& point1, const Point& point2, double unvalidBuffer/* = 0.0000001*/){ Envelope ev ; ev.expandToInclude(point1); ev.expandToInclude(point2); if(ev.getHeight() < ConstValue::ZERO || ev.getWidth() < ConstValue::ZERO) ev.expandBy(unvalidBuffer); return ev;}Envelope GeoSpatialCaculation:: MakeEnvelope(const Point& point1, double buffer){ Envelope ev; ev.expandToInclude(point1); ev.expandBy(buffer); return ev;}