LineIntersector.h

/**********************************************************************
 *
 * GEOS - Geometry Engine Open Source
 * http://geos.osgeo.org
 *
 * Copyright (C) 2005-2006 Refractions Research Inc.
 * Copyright (C) 2001-2002 Vivid Solutions Inc.
 *
 * This is free software; you can redistribute and/or modify it under
 * the terms of the GNU Lesser General Public Licence as published
 * by the Free Software Foundation.
 * See the COPYING file for more information.
 *
 **********************************************************************
 *
 * Last port: algorithm/RobustLineIntersector.java r785 (JTS-1.13+)
 *
 **********************************************************************/

#pragma once

#include <geos/export.h>
#include <geos/algorithm/Intersection.h>
#include <geos/algorithm/Interpolate.h>
#include <geos/algorithm/Orientation.h>
#include <geos/geom/Coordinate.h>
#include <geos/geom/Envelope.h>
#include <geos/geom/PrecisionModel.h>

#include <string>

// Forward declarations
namespace geos {
namespace geom {
class PrecisionModel;
}
}

namespace geos {
namespace algorithm { // geos::algorithm

class GEOS_DLL LineIntersector {
public:

    static double computeEdgeDistance(const geom::CoordinateXY& p, const geom::CoordinateXY& p0, const geom::CoordinateXY& p1);

    static double nonRobustComputeEdgeDistance(const geom::Coordinate& p, const geom::Coordinate& p1,
            const geom::Coordinate& p2);

    explicit LineIntersector(const geom::PrecisionModel* initialPrecisionModel = nullptr)
        :
        precisionModel(initialPrecisionModel),
        result(0),
        inputLines(),
        isProperVar(false)
    {}

    ~LineIntersector() = default;

    bool isInteriorIntersection()
    {
        if(isInteriorIntersection(0)) {
            return true;
        }
        if(isInteriorIntersection(1)) {
            return true;
        }
        return false;
    };

    bool isInteriorIntersection(std::size_t inputLineIndex)
    {
        for(std::size_t i = 0; i < result; ++i) {
            if(!(intPt[i].equals2D(*inputLines[inputLineIndex][0])
                 || intPt[i].equals2D(*inputLines[inputLineIndex][1]))) {
                return true;
            }
        }
        return false;
    };

    void
    setPrecisionModel(const geom::PrecisionModel* newPM)
    {
        precisionModel = newPM;
    }

    enum intersection_type : uint8_t {
        NO_INTERSECTION = 0,

        POINT_INTERSECTION = 1,

        COLLINEAR_INTERSECTION = 2
    };

    template<typename C1, typename C2>
    void computeIntersection(const C1& p1, const C1& p2,
                             const C2& p3, const C2& p4)
    {
        inputLines[0][0] = &p1;
        inputLines[0][1] = &p2;
        inputLines[1][0] = &p3;
        inputLines[1][1] = &p4;
        result = computeIntersect(p1, p2, p3, p4);
    }

    void computeIntersection(const geom::CoordinateSequence& p, std::size_t p0,
                             const geom::CoordinateSequence& q, std::size_t q0);

    std::string toString() const;

    bool
    hasIntersection() const
    {
        return result != NO_INTERSECTION;
    }


    const geom::CoordinateXY*
    getEndpoint(std::size_t segmentIndex, std::size_t ptIndex) const
    {
        return inputLines[segmentIndex][ptIndex];
    }

    size_t
    getIntersectionNum() const
    {
        return result;
    }


    const geom::CoordinateXYZM&
    getIntersection(std::size_t intIndex) const
    {
        return intPt[intIndex];
    }

    static bool isSameSignAndNonZero(double a, double b);

    bool isIntersection(const geom::Coordinate& pt) const
    {
        for(std::size_t i = 0; i < result; ++i) {
            if(intPt[i].equals2D(pt)) {
                return true;
            }
        }
        return false;
    };

    bool
    isProper() const
    {
        return hasIntersection() && isProperVar;
    }

    const geom::Coordinate& getIntersectionAlongSegment(std::size_t segmentIndex, std::size_t intIndex);

    std::size_t getIndexAlongSegment(std::size_t segmentIndex, std::size_t intIndex);

    double getEdgeDistance(std::size_t geomIndex, std::size_t intIndex) const;

private:

    const geom::PrecisionModel* precisionModel;

    std::size_t result;

    const geom::CoordinateXY* inputLines[2][2];

    geom::CoordinateXYZM intPt[2];

    std::size_t intLineIndex[2][2];

    bool isProperVar;
    //Coordinate &pa;
    //Coordinate &pb;

    bool
    isCollinear() const
    {
        return result == COLLINEAR_INTERSECTION;
    }

    template<typename C1, typename C2>
    uint8_t computeIntersect(const C1& p1, const C1& p2, const C2& q1, const C2& q2)
    {
        isProperVar = false;

        // first try a fast test to see if the envelopes of the lines intersect
        if(!geom::Envelope::intersects(p1, p2, q1, q2)) {
            return NO_INTERSECTION;
        }

        // for each endpoint, compute which side of the other segment it lies
        // if both endpoints lie on the same side of the other segment,
        // the segments do not intersect
        int Pq1 = Orientation::index(p1, p2, q1);
        int Pq2 = Orientation::index(p1, p2, q2);

        if((Pq1 > 0 && Pq2 > 0) || (Pq1 < 0 && Pq2 < 0)) {
            return NO_INTERSECTION;
        }

        int Qp1 = Orientation::index(q1, q2, p1);
        int Qp2 = Orientation::index(q1, q2, p2);

        if((Qp1 > 0 && Qp2 > 0) || (Qp1 < 0 && Qp2 < 0)) {
            return NO_INTERSECTION;
        }

        bool collinear = Pq1 == 0 && Pq2 == 0 && Qp1 == 0 && Qp2 == 0;
        if(collinear) {
            return computeCollinearIntersection(p1, p2, q1, q2);
        }

        /*
         * At this point we know that there is a single intersection point
         * (since the lines are not collinear).
         */

        /*
         * Check if the intersection is an endpoint.
         * If it is, copy the endpoint as
         * the intersection point. Copying the point rather than
         * computing it ensures the point has the exact value,
         * which is important for robustness. It is sufficient to
         * simply check for an endpoint which is on the other line,
         * since at this point we know that the inputLines must
         *  intersect.
         */
        geom::CoordinateXYZM p;
        double z = DoubleNotANumber;
        double m = DoubleNotANumber;

        if(Pq1 == 0 || Pq2 == 0 || Qp1 == 0 || Qp2 == 0) {

        isProperVar = false;

        /* Check for two equal endpoints.
         * This is done explicitly rather than by the orientation tests
         * below in order to improve robustness.
         *
         * (A example where the orientation tests fail
         *  to be consistent is:
         *
         * LINESTRING ( 19.850257749638203 46.29709338043669,
         *                      20.31970698357233 46.76654261437082 )
         * and
         * LINESTRING ( -48.51001596420236 -22.063180333403878,
         *                      19.850257749638203 46.29709338043669 )
         *
         * which used to produce the INCORRECT result:
         * (20.31970698357233, 46.76654261437082, NaN)
         */

        if (p1.equals2D(q1)) {
            p = p1;
            z = Interpolate::zGet(p1, q1);
            m = Interpolate::mGet(p1, q1);
        }
        else if (p1.equals2D(q2)) {
            p = p1;
            z = Interpolate::zGet(p1, q2);
            m = Interpolate::mGet(p1, q2);
        }
        else if (p2.equals2D(q1)) {
            p = p2;
            z = Interpolate::zGet(p2, q1);
            m = Interpolate::mGet(p2, q1);
        }
        else if (p2.equals2D(q2)) {
            p = p2;
            z = Interpolate::zGet(p2, q2);
            m = Interpolate::mGet(p2, q2);
        }
        /*
         * Now check to see if any endpoint lies on the interior of the other segment.
         */
        else if(Pq1 == 0) {
            p = q1;
            z = Interpolate::zGetOrInterpolate(q1, p1, p2);
            m = Interpolate::mGetOrInterpolate(q1, p1, p2);
        }
        else if(Pq2 == 0) {
            p = q2;
            z = Interpolate::zGetOrInterpolate(q2, p1, p2);
            m = Interpolate::mGetOrInterpolate(q2, p1, p2);
        }
        else if(Qp1 == 0) {
            p = p1;
            z = Interpolate::zGetOrInterpolate(p1, q1, q2);
            m = Interpolate::mGetOrInterpolate(p1, q1, q2);
        }
        else if(Qp2 == 0) {
            p = p2;
            z = Interpolate::zGetOrInterpolate(p2, q1, q2);
            m = Interpolate::mGetOrInterpolate(p2, q1, q2);
        }
    } else {
        isProperVar = true;
        p = intersection(p1, p2, q1, q2);
        z = Interpolate::zInterpolate(p, p1, p2, q1, q2);
        m = Interpolate::mInterpolate(p, p1, p2, q1, q2);
        }
        intPt[0] = geom::CoordinateXYZM(p.x, p.y, z, m);
    #if GEOS_DEBUG
        std::cerr << " POINT_INTERSECTION; intPt[0]:" << intPt[0].toString() << std::endl;
    #endif // GEOS_DEBUG
        return POINT_INTERSECTION;
    }

    bool
    isEndPoint() const
    {
        return hasIntersection() && !isProperVar;
    }

    void computeIntLineIndex();

    void computeIntLineIndex(std::size_t segmentIndex);

    template<typename C1, typename C2>
    uint8_t computeCollinearIntersection(const C1& p1, const C1& p2, const C2& q1, const C2& q2)
    {
        bool q1inP = geom::Envelope::intersects(p1, p2, q1);
        bool q2inP = geom::Envelope::intersects(p1, p2, q2);
        bool p1inQ = geom::Envelope::intersects(q1, q2, p1);
        bool p2inQ = geom::Envelope::intersects(q1, q2, p2);

        if(q1inP && q2inP) {
            intPt[0] = zmGetOrInterpolateCopy(q1, p1, p2);
            intPt[1] = zmGetOrInterpolateCopy(q2, p1, p2);
            return COLLINEAR_INTERSECTION;
        }
        if(p1inQ && p2inQ) {
            intPt[0] = zmGetOrInterpolateCopy(p1, q1, q2);
            intPt[1] = zmGetOrInterpolateCopy(p2, q1, q2);
            return COLLINEAR_INTERSECTION;
        }
        if(q1inP && p1inQ) {
            // if pts are equal Z is chosen arbitrarily
            intPt[0] = zmGetOrInterpolateCopy(q1, p1, p2);
            intPt[1] = zmGetOrInterpolateCopy(p1, q1, q2);

            return (q1 == p1) && !q2inP && !p2inQ ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
        }
        if(q1inP && p2inQ) {
            // if pts are equal Z is chosen arbitrarily
            intPt[0] = zmGetOrInterpolateCopy(q1, p1, p2);
            intPt[1] = zmGetOrInterpolateCopy(p2, q1, q2);

            return (q1 == p2) && !q2inP && !p1inQ ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
        }
        if(q2inP && p1inQ) {
            // if pts are equal Z is chosen arbitrarily
            intPt[0] = zmGetOrInterpolateCopy(q2, p1, p2);
            intPt[1] = zmGetOrInterpolateCopy(p1, q1, q2);

            return (q2 == p1) && !q1inP && !p2inQ ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
        }
        if(q2inP && p2inQ) {
            // if pts are equal Z is chosen arbitrarily
            intPt[0] = zmGetOrInterpolateCopy(q2, p1, p2);
            intPt[1] = zmGetOrInterpolateCopy(p2, q1, q2);
            return (q2 == p2) && !q1inP && !p1inQ ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
        }
        return NO_INTERSECTION;
    }

    template<typename C1, typename C2>
    geom::CoordinateXYZM intersection (const C1& p1, const C1& p2, const C2& q1, const C2& q2) const {
        auto intPtOut = intersectionSafe(p1, p2, q1, q2);

        /*
         * Due to rounding it can happen that the computed intersection is
         * outside the envelopes of the input segments.  Clearly this
         * is inconsistent.
         * This code checks this condition and forces a more reasonable answer
         *
         * MD - May 4 2005 - This is still a problem.  Here is a failure case:
         *
         * LINESTRING (2089426.5233462777 1180182.3877339689,
         *             2085646.6891757075 1195618.7333999649)
         * LINESTRING (1889281.8148903656 1997547.0560044837,
         *             2259977.3672235999 483675.17050843034)
         * int point = (2097408.2633752143,1144595.8008114607)
         */

        if(! isInSegmentEnvelopes(intPtOut)) {
            //intPt = CentralEndpointIntersector::getIntersection(p1, p2, q1, q2);
            intPtOut = nearestEndpoint(p1, p2, q1, q2);
        }

        if(precisionModel != nullptr) {
            precisionModel->makePrecise(intPtOut);
        }

        return intPtOut;
    }

    bool isInSegmentEnvelopes(const geom::CoordinateXY& pt) const
    {
        geom::Envelope env0(*inputLines[0][0], *inputLines[0][1]);
        geom::Envelope env1(*inputLines[1][0], *inputLines[1][1]);
        return env0.contains(pt) && env1.contains(pt);
    };

    template<typename C1, typename C2>
    geom::CoordinateXYZM intersectionSafe(const C1& p1, const C1& p2,
                                          const C2& q1, const C2& q2) const
    {
        geom::CoordinateXYZM ptInt(Intersection::intersection(p1, p2, q1, q2));
        if (ptInt.isNull()) {
            const geom::CoordinateXY& nearest = nearestEndpoint(p1, p2, q1, q2);
#if __cplusplus >= 201703L
            if constexpr (std::is_same<C1, C2>::value) {
#else
            if (std::is_same<C1, C2>::value) {
#endif
                ptInt = static_cast<const C1&>(nearest);
            } else {
                if (&nearest == static_cast<const geom::CoordinateXY*>(&p1) || &nearest == static_cast<const geom::CoordinateXY*>(&p2)) {
                    ptInt = static_cast<const C1&>(nearest);
                } else {
                    ptInt = static_cast<const C2&>(nearest);
                }
            }
        }
        return ptInt;
    }

    static const geom::CoordinateXY& nearestEndpoint(const geom::CoordinateXY& p1,
                                                     const geom::CoordinateXY& p2,
                                                     const geom::CoordinateXY& q1,
                                                     const geom::CoordinateXY& q2);


    template<typename C1, typename C2>
    static geom::CoordinateXYZM zmGetOrInterpolateCopy(
        const C1& p,
        const C2& p1,
        const C2& p2)
    {
        geom::CoordinateXYZM pCopy(p);
        pCopy.z = Interpolate::zGetOrInterpolate(p, p1, p2);
        pCopy.m = Interpolate::mGetOrInterpolate(p, p1, p2);
        return pCopy;
    }

};


} // namespace geos::algorithm
} // namespace geos