/src/geos/include/geos/algorithm/LineIntersector.h

Source
/**********************************************************************
 *
 * 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

/** \brief
 * A LineIntersector is an algorithm that can both test whether
 * two line segments intersect and compute the intersection point
 * if they do.
 *
 * The intersection point may be computed in a precise or non-precise manner.
 * Computing it precisely involves rounding it to an integer.  (This assumes
 * that the input coordinates have been made precise by scaling them to
 * an integer grid.)
 *
 */
class GEOS_DLL LineIntersector {
public:

    /// Computes the "edge distance" of an intersection point p in an edge.
    ///
    /// The edge distance is a metric of the point along the edge.
    /// The metric used is a robust and easy to compute metric function.
    /// It is <b>not</b> equivalent to the usual Euclidean metric.
    /// It relies on the fact that either the x or the y ordinates of the
    /// points in the edge are unique, depending on whether the edge is longer in
    /// the horizontal or vertical direction.
    ///
    /// NOTE: This function may produce incorrect distances
    ///  for inputs where p is not precisely on p1-p2
    /// (E.g. p = (139,9) p1 = (139,10), p2 = (280,1) produces distanct
    /// 0.0, which is incorrect.
    ///
    /// My hypothesis is that the function is safe to use for points which are the
    /// result of <b>rounding</b> points which lie on the line,
    /// but not safe to use for <b>truncated</b> points.
    ///
    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;

    /** \brief
     * Tests whether either intersection point is an interior point of
     * one of the input segments.
     *
     * @return <code>true</code> if either intersection point is in
     * the interior of one of the input segments
     */
    bool isInteriorIntersection()
    {
        if(isInteriorIntersection(0)) {
            return true;
        }
        if(isInteriorIntersection(1)) {
            return true;
        }
        return false;
    };

    /** \brief
     * Tests whether either intersection point is an interior point
     * of the specified input segment.
     *
     * @return <code>true</code> if either intersection point is in
     * the interior of the input segment
     */
    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;
    };

    /// Force computed intersection to be rounded to a given precision model.
    ///
    /// No getter is provided, because the precision model is not required
    /// to be specified.
    /// @param newPM the PrecisionModel to use for rounding
    ///
    void
    setPrecisionModel(const geom::PrecisionModel* newPM)
    {
        precisionModel = newPM;
    }

    enum intersection_type : uint8_t {
        /// Indicates that line segments do not intersect
        NO_INTERSECTION = 0,

        /// Indicates that line segments intersect in a single point
        POINT_INTERSECTION = 1,

        /// Indicates that line segments intersect in a line segment
        COLLINEAR_INTERSECTION = 2
    };

    /// Computes the intersection of the lines p1-p2 and p3-p4
    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);
    }

    /// Compute the intersection between two segments, given a sequence and starting index of each
    void computeIntersection(const geom::CoordinateSequence& p, std::size_t p0,
                             const geom::CoordinateSequence& q, std::size_t q0);

    std::string toString() const;

    /**
     * Tests whether the input geometries intersect.
     *
     * @return true if the input geometries intersect
     */
    bool
    hasIntersection() const
    {
        return result != NO_INTERSECTION;
    }


    /**
    * Gets an endpoint of an input segment.
    *
    * @param segmentIndex the index of the input segment (0 or 1)
    * @param ptIndex the index of the endpoint (0 or 1)
    * @return the specified endpoint
    */
    const geom::CoordinateXY*
    getEndpoint(std::size_t segmentIndex, std::size_t ptIndex) const
    {
        return inputLines[segmentIndex][ptIndex];
    }

    /// Returns the number of intersection points found.
    ///
    /// This will be either 0, 1 or 2.
    ///
    size_t
    getIntersectionNum() const
    {
        return result;
    }


    /// Returns the intIndex'th intersection point
    ///
    /// @param intIndex is 0 or 1
    ///
    /// @return the intIndex'th intersection point
    ///
    const geom::CoordinateXYZM&
    getIntersection(std::size_t intIndex) const
    {
        return intPt[intIndex];
    }

    /// Returns false if both numbers are zero.
    ///
    /// @return true if both numbers are positive or if both numbers are negative.
    ///
    static bool isSameSignAndNonZero(double a, double b);

    /** \brief
     * Test whether a point is a intersection point of two line segments.
     *
     * Note that if the intersection is a line segment, this method only tests for
     * equality with the endpoints of the intersection segment.
     * It does <b>not</b> return true if
     * the input point is internal to the intersection segment.
     *
     * @return true if the input point is one of the intersection points.
     */
    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;
    };

    /** \brief
     * Tests whether an intersection is proper.
     *
     * The intersection between two line segments is considered proper if
     * they intersect in a single point in the interior of both segments
     * (e.g. the intersection is a single point and is not equal to any of the
     * endpoints).
     *
     * The intersection between a point and a line segment is considered proper
     * if the point lies in the interior of the segment (e.g. is not equal to
     * either of the endpoints).
     *
     * @return true if the intersection is proper
     */
    bool
    isProper() const
    {
        return hasIntersection() && isProperVar;
    }

    /** \brief
     * Computes the intIndex'th intersection point in the direction of
     * a specified input line segment
     *
     * @param segmentIndex is 0 or 1
     * @param intIndex is 0 or 1
     *
     * @return the intIndex'th intersection point in the direction of the
     *         specified input line segment
     */
    const geom::Coordinate& getIntersectionAlongSegment(std::size_t segmentIndex, std::size_t intIndex);

    /** \brief
     * Computes the index of the intIndex'th intersection point in the direction of
     * a specified input line segment
     *
     * @param segmentIndex is 0 or 1
     * @param intIndex is 0 or 1
     *
     * @return the index of the intersection point along the segment (0 or 1)
     */
    std::size_t getIndexAlongSegment(std::size_t segmentIndex, std::size_t intIndex);

    /** \brief
     * Computes the "edge distance" of an intersection point along the specified
     * input line segment.
     *
     * @param geomIndex is 0 or 1
     * @param intIndex is 0 or 1
     *
     * @return the edge distance of the intersection point
     */
    double getEdgeDistance(std::size_t geomIndex, std::size_t intIndex) const;

private:

    /**
     * If makePrecise is true, computed intersection coordinates
     * will be made precise using Coordinate#makePrecise
     */
    const geom::PrecisionModel* precisionModel;

    std::size_t result;

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

    /**
     * We store real Coordinates here because
     * we must compute the Z of intersection point.
     */
    geom::CoordinateXYZM intPt[2];

    /**
     * The indexes of the endpoints of the intersection lines, in order along
     * the corresponding line
     */
    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;
        }

        /**
         * Intersection is collinear if each endpoint lies on the other line.
         */
        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;
    }

    /** \brief
     * This method computes the actual value of the intersection point.
     *
     * To obtain the maximum precision from the intersection calculation,
     * the coordinates are normalized by subtracting the minimum
     * ordinate values (in absolute value).  This has the effect of
     * removing common significant digits from the calculation to
     * maintain more bits of precision.
     */
    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;
    }

    /**
     * Test whether a point lies in the envelopes of both input segments.
     * A correctly computed intersection point should return true
     * for this test.
     * Since this test is for debugging purposes only, no attempt is
     * made to optimize the envelope test.
     *
     * @return true if the input point lies within both
     *         input segment envelopes
     */
    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);
    };

    /**
     * Computes a segment intersection.
     * Round-off error can cause the raw computation to fail,
     * (usually due to the segments being approximately parallel).
     * If this happens, a reasonable approximation is computed instead.
     *
     * @param p1 a segment endpoint
     * @param p2 a segment endpoint
     * @param q1 a segment endpoint
     * @param q2 a segment endpoint
     * @return the computed intersection point is stored there
     */
    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;
    }

    /**
     * Finds the endpoint of the segments P and Q which
     * is closest to the other segment.
     * This is a reasonable surrogate for the true
     * intersection points in ill-conditioned cases
     * (e.g. where two segments are nearly coincident,
     * or where the endpoint of one segment lies almost on the other segment).
     * <p>
     * This replaces the older CentralEndpoint heuristic,
     * which chose the wrong endpoint in some cases
     * where the segments had very distinct slopes
     * and one endpoint lay almost on the other segment.
     *
     * @param p1 an endpoint of segment P
     * @param p2 an endpoint of segment P
     * @param q1 an endpoint of segment Q
     * @param q2 an endpoint of segment Q
     * @return the nearest endpoint to the other segment
     */
    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

















Coverage Report

Created: 2026-01-22 06:39