From 5cbfdb39240296dbe2f9858c52521e35b3c43fe5 Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Fri, 17 Nov 2023 10:23:45 -0800 Subject: [PATCH] Port JTS-819, Add IntersectionLineSegment function, Closes GH-995, Buffer with mitre join produces error: orientationIndex encountered NaN/Inf numbers --- include/geos/algorithm/Intersection.h | 51 ++++++-- src/algorithm/Intersection.cpp | 63 ++++++++- .../buffer/OffsetSegmentGenerator.cpp | 29 +---- tests/unit/algorithm/IntersectionTest.cpp | 120 +++++++++++++++--- tests/unit/capi/GEOSBufferTest.cpp | 15 +++ 5 files changed, 218 insertions(+), 60 deletions(-) diff --git a/include/geos/algorithm/Intersection.h b/include/geos/algorithm/Intersection.h index c4a38fb4f0..de413d4a10 100644 --- a/include/geos/algorithm/Intersection.h +++ b/include/geos/algorithm/Intersection.h @@ -19,19 +19,28 @@ #include namespace geos { -namespace algorithm { // geos::algorithm +namespace algorithm { + + +/** \brief + * Functions to compute intersection points between lines and line segments. + * + * In general it is not possible to compute + * the intersection point of two lines exactly, due to numerical roundoff. + * This is particularly true when the lines are nearly parallel. + * These routines uses numerical conditioning on the input values + * to ensure that the computed value is very close to the correct value. + * + * The Z-ordinate is ignored, and not populated. + */ +class GEOS_DLL Intersection { + +public: /** \brief * Computes the intersection point of two lines. * If the lines are parallel or collinear this case is detected * and null is returned. - *

- * In general it is not possible to accurately compute - * the intersection point of two lines, due to - * numerical roundoff. - * This is particularly true when the input lines are nearly parallel. - * This routine uses numerical conditioning on the input values - * to ensure that the computed value should be very close to the correct value. * * @param p1 an endpoint of line 1 * @param p2 an endpoint of line 1 @@ -42,13 +51,31 @@ namespace algorithm { // geos::algorithm * * @see CGAlgorithmsDD#intersection(Coordinate, Coordinate, Coordinate, Coordinate) */ -class GEOS_DLL Intersection { - -public: - static geom::CoordinateXY intersection(const geom::CoordinateXY& p1, const geom::CoordinateXY& p2, const geom::CoordinateXY& q1, const geom::CoordinateXY& q2); +/** +* Computes the intersection point of a line and a line segment (if any). +* There will be no intersection point if: +* +* * the segment does not intersect the line +* * the line or the segment are degenerate (have zero length) +* +* If the segment is collinear with the line the first segment endpoint is returned. +* +* @param line1 a point on the line +* @param line2 a point on the line +* @param seg1 an endpoint of the line segment +* @param seg2 an endpoint of the line segment +* @return the intersection point, or null if it is not possible to find an intersection +*/ +static geom::CoordinateXY intersectionLineSegment( + const geom::CoordinateXY& line1, + const geom::CoordinateXY& line2, + const geom::CoordinateXY& seg1, + const geom::CoordinateXY& seg2); + + }; diff --git a/src/algorithm/Intersection.cpp b/src/algorithm/Intersection.cpp index e23c33dccb..7970f6f8fa 100644 --- a/src/algorithm/Intersection.cpp +++ b/src/algorithm/Intersection.cpp @@ -15,8 +15,10 @@ #include #include -#include #include +#include +#include +#include namespace geos { namespace algorithm { // geos.algorithm @@ -31,6 +33,65 @@ Intersection::intersection(const geom::CoordinateXY& p1, const geom::CoordinateX } +/** +* Computes the intersection point of a line and a line segment (if any). +* There will be no intersection point if: +* +* * the segment does not intersect the line +* * the line or the segment are degenerate (have zero length) +* +* If the segment is collinear with the line the first segment endpoint is returned. +* +* @param line1 a point on the line +* @param line2 a point on the line +* @param seg1 an endpoint of the line segment +* @param seg2 an endpoint of the line segment +* @return the intersection point, or null if it is not possible to find an intersection +*/ +/* public static */ +geom::CoordinateXY +Intersection::intersectionLineSegment( + const geom::CoordinateXY& line1, + const geom::CoordinateXY& line2, + const geom::CoordinateXY& seg1, + const geom::CoordinateXY& seg2) +{ + int orientS1 = Orientation::index(line1, line2, seg1); + if (orientS1 == 0) + return seg1; + + int orientS2 = Orientation::index(line1, line2, seg2); + if (orientS2 == 0) + return seg2; + + /** + * If segment lies completely on one side of the line, it does not intersect + */ + if ((orientS1 > 0 && orientS2 > 0) || (orientS1 < 0 && orientS2 < 0)) { + return geom::CoordinateXY::getNull(); + } + + /** + * The segment intersects the line. + * The full line-line intersection is used to compute the intersection point. + */ + geom::CoordinateXY intPt = intersection(line1, line2, seg1, seg2); + if (!intPt.isNull()) + return intPt; + + /** + * Due to robustness failure it is possible the intersection computation will return null. + * In this case choose the closest point + */ + double dist1 = Distance::pointToLinePerpendicular(seg1, line1, line2); + double dist2 = Distance::pointToLinePerpendicular(seg2, line1, line2); + if (dist1 < dist2) + return seg1; + else + return seg2; +} + + } // namespace geos.algorithm } // namespace geos diff --git a/src/operation/buffer/OffsetSegmentGenerator.cpp b/src/operation/buffer/OffsetSegmentGenerator.cpp index 915bbd9548..4537f9b3d5 100644 --- a/src/operation/buffer/OffsetSegmentGenerator.cpp +++ b/src/operation/buffer/OffsetSegmentGenerator.cpp @@ -35,10 +35,6 @@ #include #include -#ifndef GEOS_DEBUG -#define GEOS_DEBUG 0 -#endif - using namespace geos::algorithm; using namespace geos::geom; @@ -538,15 +534,9 @@ OffsetSegmentGenerator::addLimitedMitreJoin( // compute the candidate bevel segment by projecting both sides of the midpoint Coordinate bevel0 = project(bevelMidPt, p_distance, dirBevel); Coordinate bevel1 = project(bevelMidPt, p_distance, dirBevel + MATH_PI); - LineSegment bevel(bevel0, bevel1); - - //-- compute intersections with extended offset segments - double extendLen = p_mitreLimitDistance < p_distance ? p_distance : p_mitreLimitDistance; - LineSegment extend0 = extend(p_offset0, 2 * extendLen); - LineSegment extend1 = extend(p_offset1, -2 * extendLen); - Coordinate bevelInt0 = bevel.intersection(extend0); - Coordinate bevelInt1 = bevel.intersection(extend1); + Coordinate bevelInt0(Intersection::intersectionLineSegment(p_offset0.p0, p_offset0.p1, bevel0, bevel1)); + Coordinate bevelInt1(Intersection::intersectionLineSegment(p_offset1.p0, p_offset1.p1, bevel0, bevel1)); //-- add the limited bevel, if it intersects the offsets if (!bevelInt0.isNull() && !bevelInt1.isNull()) { @@ -563,21 +553,6 @@ OffsetSegmentGenerator::addLimitedMitreJoin( } -/* private static */ -LineSegment -OffsetSegmentGenerator::extend(const LineSegment& seg, double dist) -{ - double distFrac = std::abs(dist) / seg.getLength(); - double segFrac = dist >= 0 ? 1 + distFrac : 0 - distFrac; - Coordinate extendPt; - seg.pointAlong(segFrac, extendPt); - if (dist > 0) - return LineSegment(seg.p0, extendPt); - else - return LineSegment(extendPt, seg.p1); -} - - /* private static */ Coordinate OffsetSegmentGenerator::project(const Coordinate& pt, double d, double dir) diff --git a/tests/unit/algorithm/IntersectionTest.cpp b/tests/unit/algorithm/IntersectionTest.cpp index 75adfce942..b2b6576e27 100644 --- a/tests/unit/algorithm/IntersectionTest.cpp +++ b/tests/unit/algorithm/IntersectionTest.cpp @@ -43,8 +43,10 @@ struct test_intersection_data { double MAX_ABS_ERROR = 1e-5; - void checkIntersectionNull(double p1x, double p1y, double p2x, double p2y, - double q1x, double q1y, double q2x, double q2y) { + void checkIntersectionNull( + double p1x, double p1y, double p2x, double p2y, + double q1x, double q1y, double q2x, double q2y) + { Coordinate p1(p1x, p1y); Coordinate p2(p2x, p2y); Coordinate q1(q1x, q1y); @@ -53,9 +55,11 @@ struct test_intersection_data { ensure("checkIntersectionNull", actual.isNull()); } - void checkIntersection(double p1x, double p1y, double p2x, double p2y, - double q1x, double q1y, double q2x, double q2y, - double expectedx, double expectedy) { + void checkIntersection( + double p1x, double p1y, double p2x, double p2y, + double q1x, double q1y, double q2x, double q2y, + double expectedx, double expectedy) + { Coordinate p1(p1x, p1y); Coordinate p2(p2x, p2y); Coordinate q1(q1x, q1y); @@ -67,6 +71,33 @@ struct test_intersection_data { ensure("checkIntersection", dist <= MAX_ABS_ERROR); } + void checkIntersectionLineSegment( + double p1x, double p1y, double p2x, double p2y, + double q1x, double q1y, double q2x, double q2y, + double expectedx, double expectedy) + { + Coordinate p1(p1x, p1y); + Coordinate p2(p2x, p2y); + Coordinate q1(q1x, q1y); + Coordinate q2(q2x, q2y); + Coordinate actual(Intersection::intersectionLineSegment(p1, p2, q1, q2)); + Coordinate expected(expectedx, expectedy); + double dist = actual.distance(expected); + ensure("checkIntersectionLineSegment", dist <= MAX_ABS_ERROR); + } + + void checkIntersectionLineSegmentNull( + double p1x, double p1y, double p2x, double p2y, + double q1x, double q1y, double q2x, double q2y) + { + Coordinate p1(p1x, p1y); + Coordinate p2(p2x, p2y); + Coordinate q1(q1x, q1y); + Coordinate q2(q2x, q2y); + Coordinate actual(Intersection::intersectionLineSegment(p1, p2, q1, q2)); + ensure("checkIntersectionLineSegmentNull", actual.isNull()); + } + test_intersection_data() {} @@ -77,41 +108,90 @@ typedef group::object object; group test_intersection_data("geos::algorithm::Intersection"); - +// testSimple template<> template<> -void object::test<1> -() +void object::test<1>() { - // testSimple - checkIntersection( - 0,0, 10,10, - 0,10, 10,0, - 5,5); + checkIntersection(0,0, 10,10, 0,10, 10,0, 5,5); +} - // testCollinear - checkIntersectionNull( - 0,0, 10,10, - 20,20, 30, 30 ); +// testCollinear +template<> +template<> +void object::test<2>() +{ + checkIntersectionNull(0,0, 10,10, 20,20, 30, 30 ); +} - // testParallel +// testParallel +template<> +template<> +void object::test<3>() +{ checkIntersectionNull( 0,0, 10,10, 10,0, 20,10 ); +} - // testAlmostCollinear +// testAlmostCollinear +template<> +template<> +void object::test<4>() +{ checkIntersection( 35613471.6165017, 4257145.306132293, 35613477.7705378, 4257160.528222711, 35613477.77505724, 4257160.539653536, 35613479.85607389, 4257165.92369170, 35613477.772841461, 4257160.5339209242 ); +} - // testAlmostCollinearCond +// testAlmostCollinearCond +template<> +template<> +void object::test<5>() +{ checkIntersection( 1.6165017, 45.306132293, 7.7705378, 60.528222711, 7.77505724, 60.539653536, 9.85607389, 65.92369170, 7.772841461, 60.5339209242 ); +} +// testLineSegCross +template<> +template<> +void object::test<6>() +{ + checkIntersectionLineSegment( 0, 0, 0, 1, -1, 9, 1, 9, 0, 9 ); + checkIntersectionLineSegment( 0, 0, 0, 1, -1, 2, 1, 4, 0, 3 ); } +// testLineSegTouch +template<> +template<> +void object::test<7>() +{ + checkIntersectionLineSegment( 0, 0, 0, 1, -1, 9, 0, 9, 0, 9 ); + checkIntersectionLineSegment( 0, 0, 0, 1, 0, 2, 1, 4, 0, 2 ); +} + +// testLineSegCollinear +template<> +template<> +void object::test<8>() +{ + checkIntersectionLineSegment( 0, 0, 0, 1, 0, 9, 0, 8, 0, 9 ); +} + +// testLineSegNone +template<> +template<> +void object::test<9>() +{ + checkIntersectionLineSegmentNull( 0, 0, 0, 1, 2, 9, 1, 9 ); + checkIntersectionLineSegmentNull( 0, 0, 0, 1, -2, 9, -1, 9 ); + checkIntersectionLineSegmentNull( 0, 0, 0, 1, 2, 9, 1, 9 ); +} + + } // namespace tut diff --git a/tests/unit/capi/GEOSBufferTest.cpp b/tests/unit/capi/GEOSBufferTest.cpp index f47a2237f9..e1124c2c27 100644 --- a/tests/unit/capi/GEOSBufferTest.cpp +++ b/tests/unit/capi/GEOSBufferTest.cpp @@ -615,4 +615,19 @@ void object::test<24> ensure(result_ == nullptr); } + +template<> +template<> +void object::test<25> +() +{ + geom1_ = GEOSGeomFromWKT("POLYGON ((4.6664239253667485 4.9470840685113275, 4.666423925366749 4.947084068511328, 3.569508914897422 -10.739531408188364, -9.082056557097435 19.893317266250286, 5.639581102785941 18.86388007810711, 4.6664239253667485 4.9470840685113275))"); + geom2_ = GEOSBufferWithStyle(geom1_, -1, 8, + GEOSBUF_CAP_ROUND, + GEOSBUF_JOIN_MITRE, + 5); + geom3_ = GEOSGeomFromWKT("POLYGON ((3.3225774291798533 0.0647708524944821, 3.3225774291798555 0.0647708524944812, 2.8688758567150883 -6.4234639154696263, -7.5416226086581215 18.7831577331451953, 4.5722605787819921 17.9360725015914078, 3.3225774291798533 0.0647708524944821))"); + ensure_geometry_equals_exact(geom3_, geom2_, 0.001); +} + } // namespace tut