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

@@ -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);
+
+
 };
 
 

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
 

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)

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
 

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