二维空间计算几何

Posted 2021-02-09 dingshengtao

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了二维空间计算几何相关的知识，希望对你有一定的参考价值。

分享一下之前做项目写的空间计算几何源码，

第一部分，涵盖简单计算几何方法的实现。例如：外包框计算、距离计算、面积与长度计算、空间关系的判断、方位角、交点计算、多边形中心点、多边形最大内接圆、凸包、投影点等。

第二部分主要是复杂的空间计算方法：多边形的交差并补计算、缓冲区、求平行线、求延长线等

第三部分主要是分享一下空间拓扑。

第二、三部分以后补上。

一. 头文件

//************************************** 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);

View Code

二.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<Point>& 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<Point>& points)
    {
        std::vector<Point> 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<Point>& 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<Point>& 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<Point>& pointVecIntersect, GeoVectorT<Point>& 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 
", j);
                GeoVectorT<Point> 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);
    }

View Code

//************************************** 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<int> stack;
    stack.push_back(0); 
    stack.push_back(1);

    int top=2;
    int n = pointsRet.size();
    for (int i=2;i<n;i++)
    {
        // 由S的栈顶元素的下一个元素、S的栈顶元素以及pi构成的折线段拐向右侧，S出栈
        while (top>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<stack.size();i++)//输出栈中的所有点，即为凸包上的点
    {
        if(stack[i] == 0 && i != 0)
            break;
        pointVecRet.push_back(pointsRet[stack[i]]);
    }
}
Polygon* GeoSpatialCaculation::             PointsConvexHull(const PointVec& points)
{
    PointVec pointVecRet;
    PointsConvexHull(points, pointVecRet);
    if(pointVecRet[0] != pointVecRet[pointVecRet.size()-1]) pointVecRet.push_back(pointVecRet[0]);    
    return pointVecRet.size() < 4 ? 0 : GeometryFactory::createPolygon(GeometryFactory::createLinearRing(pointVecRet), true);
}
Point GeoSpatialCaculation::                PointsCentroid(const PointVec& points)
{
    double x = 0, y = 0;
    for(int i = 0; i < points.size(); i++)
    {
        x += points[i].getX();
        y += points[i].getY();
    }
    return Point(x/points.size(), y/points.size());
}
Point GeoSpatialCaculation::                PolygonInnerPoint(const Polygon* polygon)
{
    Point pt(ConstValue::MAXVALUE, ConstValue::MAXVALUE);
    const LinearRing* linearRing = polygon->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;
}

View Code

以上是关于二维空间计算几何的主要内容，如果未能解决你的问题，请参考以下文章