[geos-commits] [SCM] GEOS branch main updated. 5cbfdb39240296dbe2f9858c52521e35b3c43fe5
git at osgeo.org
git at osgeo.org
Fri Nov 17 11:54:39 PST 2023
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "GEOS".
The branch, main has been updated
via 5cbfdb39240296dbe2f9858c52521e35b3c43fe5 (commit)
from dcde8ad8a15eabdafd3a7c3ef78d6cf20cf800de (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 5cbfdb39240296dbe2f9858c52521e35b3c43fe5
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Fri Nov 17 10:23:45 2023 -0800
Port JTS-819, Add IntersectionLineSegment function,
Closes GH-995, Buffer with mitre join produces error: orientationIndex encountered NaN/Inf numbers
diff --git a/include/geos/algorithm/Intersection.h b/include/geos/algorithm/Intersection.h
index c4a38fb4f..de413d4a1 100644
--- a/include/geos/algorithm/Intersection.h
+++ b/include/geos/algorithm/Intersection.h
@@ -19,19 +19,28 @@
#include <geos/geom/CoordinateSequence.h>
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 <code>null</code> is returned.
- * <p>
- * 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 e23c33dcc..7970f6f8f 100644
--- a/src/algorithm/Intersection.cpp
+++ b/src/algorithm/Intersection.cpp
@@ -15,8 +15,10 @@
#include <cmath>
#include <vector>
-#include <geos/algorithm/Intersection.h>
#include <geos/algorithm/CGAlgorithmsDD.h>
+#include <geos/algorithm/Distance.h>
+#include <geos/algorithm/Intersection.h>
+#include <geos/algorithm/Orientation.h>
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 915bbd954..4537f9b3d 100644
--- a/src/operation/buffer/OffsetSegmentGenerator.cpp
+++ b/src/operation/buffer/OffsetSegmentGenerator.cpp
@@ -35,10 +35,6 @@
#include <geos/algorithm/CGAlgorithmsDD.h>
#include <geos/util.h>
-#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 75adfce94..b2b6576e2 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 f47a2237f..e1124c2c2 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
-----------------------------------------------------------------------
Summary of changes:
include/geos/algorithm/Intersection.h | 51 +++++++---
src/algorithm/Intersection.cpp | 63 ++++++++++++-
src/operation/buffer/OffsetSegmentGenerator.cpp | 29 +-----
tests/unit/algorithm/IntersectionTest.cpp | 120 ++++++++++++++++++++----
tests/unit/capi/GEOSBufferTest.cpp | 15 +++
5 files changed, 218 insertions(+), 60 deletions(-)
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list