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
geos::algorithm::LineIntersector
A LineIntersector is an algorithm that can both test whether two line segments intersect and compute ...
Definition LineIntersector.h:53
geos::algorithm::LineIntersector::intersection_type
intersection_type
Definition LineIntersector.h:137
geos::algorithm::LineIntersector::NO_INTERSECTION
@ NO_INTERSECTION
Indicates that line segments do not intersect.
Definition LineIntersector.h:139
geos::algorithm::LineIntersector::COLLINEAR_INTERSECTION
@ COLLINEAR_INTERSECTION
Indicates that line segments intersect in a line segment.
Definition LineIntersector.h:145
geos::algorithm::LineIntersector::POINT_INTERSECTION
@ POINT_INTERSECTION
Indicates that line segments intersect in a single point.
Definition LineIntersector.h:142
geos::algorithm::LineIntersector::computeIntersection
void computeIntersection(const geom::CoordinateSequence &p, std::size_t p0, const geom::CoordinateSequence &q, std::size_t q0)
Compute the intersection between two segments, given a sequence and starting index of each.
geos::algorithm::LineIntersector::isInteriorIntersection
bool isInteriorIntersection(std::size_t inputLineIndex)
Tests whether either intersection point is an interior point of the specified input segment.
Definition LineIntersector.h:114
geos::algorithm::LineIntersector::getIntersectionAlongSegment
const geom::Coordinate & getIntersectionAlongSegment(std::size_t segmentIndex, std::size_t intIndex)
Computes the intIndex'th intersection point in the direction of a specified input line segment.
geos::algorithm::LineIntersector::setPrecisionModel
void setPrecisionModel(const geom::PrecisionModel *newPM)
Definition LineIntersector.h:132
geos::algorithm::LineIntersector::isInteriorIntersection
bool isInteriorIntersection()
Tests whether either intersection point is an interior point of one of the input segments.
Definition LineIntersector.h:96
geos::algorithm::LineIntersector::getIntersection
const geom::CoordinateXYZM & getIntersection(std::size_t intIndex) const
Definition LineIntersector.h:209
geos::algorithm::LineIntersector::computeEdgeDistance
static double computeEdgeDistance(const geom::CoordinateXY &p, const geom::CoordinateXY &p0, const geom::CoordinateXY &p1)
geos::algorithm::LineIntersector::hasIntersection
bool hasIntersection() const
Definition LineIntersector.h:172
geos::algorithm::LineIntersector::isSameSignAndNonZero
static bool isSameSignAndNonZero(double a, double b)
geos::algorithm::LineIntersector::computeIntersection
void computeIntersection(const C1 &p1, const C1 &p2, const C2 &p3, const C2 &p4)
Computes the intersection of the lines p1-p2 and p3-p4.
Definition LineIntersector.h:150
geos::algorithm::LineIntersector::getEndpoint
const geom::CoordinateXY * getEndpoint(std::size_t segmentIndex, std::size_t ptIndex) const
Definition LineIntersector.h:186
geos::algorithm::LineIntersector::isProper
bool isProper() const
Tests whether an intersection is proper.
Definition LineIntersector.h:255
geos::algorithm::LineIntersector::getIntersectionNum
size_t getIntersectionNum() const
Definition LineIntersector.h:196
geos::algorithm::LineIntersector::getEdgeDistance
double getEdgeDistance(std::size_t geomIndex, std::size_t intIndex) const
Computes the "edge distance" of an intersection point along the specified input line segment.
geos::algorithm::LineIntersector::getIndexAlongSegment
std::size_t getIndexAlongSegment(std::size_t segmentIndex, std::size_t intIndex)
Computes the index of the intIndex'th intersection point in the direction of a specified input line s...
geos::algorithm::LineIntersector::isIntersection
bool isIntersection(const geom::Coordinate &pt) const
Test whether a point is a intersection point of two line segments.
Definition LineIntersector.h:230
geos::geom::CoordinateSequence
The internal representation of a list of coordinates inside a Geometry.
Definition CoordinateSequence.h:56
geos::geom::Coordinate
Coordinate is the lightweight class used to store coordinates.
Definition Coordinate.h:217
geos::geom::Envelope::intersects
static bool intersects(const CoordinateXY &p1, const CoordinateXY &p2, const CoordinateXY &q)
Test the point q to see whether it intersects the Envelope defined by p1-p2.
geos::geom::PrecisionModel
Specifies the precision model of the Coordinate in a Geometry.
Definition PrecisionModel.h:87
geos::geom::PrecisionModel::makePrecise
double makePrecise(double val) const
Rounds a numeric value to the PrecisionModel grid.
geos::algorithm
Contains classes and interfaces implementing fundamental computational geometry algorithms.
Definition Angle.h:32
geos::geom
Definition Angle.h:26
geos
Basic namespace for all GEOS functionalities.
Definition geos.h:38