[geos-commits] [SCM] GEOS branch master updated. 32af81367a06cc47725cbf7d0c3b6773862d04d0
git at osgeo.org
git at osgeo.org
Fri Sep 6 13:32:38 PDT 2019
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, master has been updated
via 32af81367a06cc47725cbf7d0c3b6773862d04d0 (commit)
from 7a143371ed80b0f4a613894ab5717ca4cd878a21 (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 32af81367a06cc47725cbf7d0c3b6773862d04d0
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Fri Sep 6 13:31:50 2019 -0700
Port JTS PR 468
https://github.com/locationtech/jts/pull/468
Refactor line intersection code into Intersection class.
Use new class for all line intersection computation
Unroll input conditioning for performance
Add tests to verify accuracy of conditioned double-precision intersection algorithm
diff --git a/capi/geos_ts_c.cpp b/capi/geos_ts_c.cpp
index ad97e95..05e2678 100644
--- a/capi/geos_ts_c.cpp
+++ b/capi/geos_ts_c.cpp
@@ -6579,11 +6579,9 @@ extern "C" {
try {
geos::geom::LineSegment a(ax0, ay0, ax1, ay1);
geos::geom::LineSegment b(bx0, by0, bx1, by1);
- geos::geom::Coordinate isect;
+ geos::geom::Coordinate isect = a.intersection(b);
- bool intersects = a.intersection(b, isect);
-
- if(!intersects) {
+ if(isect.isNull()) {
return -1;
}
diff --git a/include/geos/algorithm/CGAlgorithmsDD.h b/include/geos/algorithm/CGAlgorithmsDD.h
index 258361b..f596466 100644
--- a/include/geos/algorithm/CGAlgorithmsDD.h
+++ b/include/geos/algorithm/CGAlgorithmsDD.h
@@ -111,9 +111,17 @@ public:
return CGAlgorithmsDD::STRAIGHT;
}
- static void intersection(const geom::Coordinate& p1, const geom::Coordinate& p2,
- const geom::Coordinate& q1, const geom::Coordinate& q2,
- geom::Coordinate& rv);
+ /**
+ * If the lines are parallel (either identical
+ * or separate) a null value is returned.
+ * @param p1 an endpoint of line segment 1
+ * @param p2 an endpoint of line segment 1
+ * @param q1 an endpoint of line segment 2
+ * @param q2 an endpoint of line segment 2
+ * @return an intersection point if one exists, or null if the lines are parallel
+ */
+ static geom::Coordinate intersection(const geom::Coordinate& p1, const geom::Coordinate& p2,
+ const geom::Coordinate& q1, const geom::Coordinate& q2);
static int signOfDet2x2(double dx1, double dy1, double dx2, double dy2);
diff --git a/include/geos/algorithm/Intersection.h b/include/geos/algorithm/Intersection.h
new file mode 100644
index 0000000..16991a3
--- /dev/null
+++ b/include/geos/algorithm/Intersection.h
@@ -0,0 +1,61 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2019 Paul Ramsey <pramsey at cleverlephant.ca>
+ *
+ * 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.
+ *
+ **********************************************************************/
+
+#ifndef GEOS_ALGORITHM_INTERSECTION_H
+#define GEOS_ALGORITHM_INTERSECTION_H
+
+#include <geos/export.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/CoordinateSequence.h>
+
+namespace geos {
+namespace algorithm { // geos::algorithm
+
+/**
+* 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
+* @param q1 an endpoint of line 2
+* @param q2 an endpoint of line 2
+* @return the intersection point between the lines, if there is one,
+* or null if the lines are parallel or collinear
+*
+* @see CGAlgorithmsDD#intersection(Coordinate, Coordinate, Coordinate, Coordinate)
+*/
+
+class GEOS_DLL Intersection {
+
+public:
+
+static geom::Coordinate intersection(const geom::Coordinate& p1, const geom::Coordinate& p2,
+ const geom::Coordinate& q1, const geom::Coordinate& q2);
+
+};
+
+
+} // namespace geos::algorithm
+} // namespace geos
+
+
+#endif // GEOS_ALGORITHM_INTERSECTION_H
diff --git a/include/geos/algorithm/LineIntersector.h b/include/geos/algorithm/LineIntersector.h
index fad89ad..94b041a 100644
--- a/include/geos/algorithm/LineIntersector.h
+++ b/include/geos/algorithm/LineIntersector.h
@@ -260,12 +260,6 @@ public:
private:
- void intersectionWithNormalization(const geom::Coordinate& p1,
- const geom::Coordinate& p2,
- const geom::Coordinate& q1,
- const geom::Coordinate& q2,
- geom::Coordinate& ret) const;
-
/**
* If makePrecise is true, computed intersection coordinates
* will be made precise using Coordinate#makePrecise
@@ -298,8 +292,8 @@ private:
return result == COLLINEAR_INTERSECTION;
}
- int computeIntersect(const geom::Coordinate& p1, const geom::Coordinate& p2, const geom::Coordinate& q1,
- const geom::Coordinate& q2);
+ int computeIntersect(const geom::Coordinate& p1, const geom::Coordinate& p2,
+ const geom::Coordinate& q1, const geom::Coordinate& q2);
bool
isEndPoint() const
@@ -311,9 +305,8 @@ private:
void computeIntLineIndex(int segmentIndex);
- int computeCollinearIntersection(const geom::Coordinate& p1,
- const geom::Coordinate& p2, const geom::Coordinate& q1,
- const geom::Coordinate& q2);
+ int computeCollinearIntersection(const geom::Coordinate& p1, const geom::Coordinate& p2,
+ const geom::Coordinate& q1, const geom::Coordinate& q2);
/** \brief
* This method computes the actual value of the intersection point.
@@ -324,14 +317,10 @@ private:
* removing common significant digits from the calculation to
* maintain more bits of precision.
*/
- void intersection(const geom::Coordinate& p1,
- const geom::Coordinate& p2,
- const geom::Coordinate& q1,
- const geom::Coordinate& q2,
- geom::Coordinate& ret) const;
-
- double smallestInAbsValue(double x1, double x2,
- double x3, double x4) const;
+ geom::Coordinate intersection(const geom::Coordinate& p1,
+ const geom::Coordinate& p2,
+ const geom::Coordinate& q1,
+ const geom::Coordinate& q2) const;
/**
* Test whether a point lies in the envelopes of both input segments.
@@ -345,23 +334,9 @@ private:
*/
bool isInSegmentEnvelopes(const geom::Coordinate& intPt) const;
- /** \brief
- * Normalize the supplied coordinates to
- * so that the midpoint of their intersection envelope
- * lies at the origin.
- *
- * @param n00
- * @param n01
- * @param n10
- * @param n11
- * @param normPt
- */
- void normalizeToEnvCentre(geom::Coordinate& n00, geom::Coordinate& n01,
- geom::Coordinate& n10, geom::Coordinate& n11,
- geom::Coordinate& normPt) const;
/**
- * Computes a segment intersection using homogeneous coordinates.
+ * 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.
@@ -370,13 +345,10 @@ private:
* @param p2 a segment endpoint
* @param q1 a segment endpoint
* @param q2 a segment endpoint
- * @param intPt the computed intersection point is stored there
+ * @return the computed intersection point is stored there
*/
- void safeHCoordinateIntersection(const geom::Coordinate& p1,
- const geom::Coordinate& p2,
- const geom::Coordinate& q1,
- const geom::Coordinate& q2,
- geom::Coordinate& intPt) const;
+ geom::Coordinate intersectionSafe(const geom::Coordinate& p1, const geom::Coordinate& p2,
+ const geom::Coordinate& q1, const geom::Coordinate& q2) const;
};
diff --git a/include/geos/algorithm/Makefile.am b/include/geos/algorithm/Makefile.am
index e4c22f1..717e63d 100644
--- a/include/geos/algorithm/Makefile.am
+++ b/include/geos/algorithm/Makefile.am
@@ -24,6 +24,7 @@ geos_HEADERS = \
InteriorPointArea.h \
InteriorPointLine.h \
InteriorPointPoint.h \
+ Intersection.h \
Length.h \
LineIntersector.h \
MinimumBoundingCircle.h \
diff --git a/include/geos/geom/LineSegment.h b/include/geos/geom/LineSegment.h
index 2430c12..38191d9 100644
--- a/include/geos/geom/LineSegment.h
+++ b/include/geos/geom/LineSegment.h
@@ -335,10 +335,9 @@ public:
* intersection, the LineIntersector class should be used.
*
* @param line
- * @param coord the Coordinate to write the result into
- * @return true if an intersection was found, false otherwise
+ * @return intersection if found, setNull() otherwise
*/
- bool intersection(const LineSegment& line, Coordinate& coord) const;
+ Coordinate intersection(const LineSegment& line) const;
/** \brief
* Computes the intersection point of the lines defined
@@ -353,11 +352,10 @@ public:
* be used.
*
* @param line a line segment defining a straight line
- * @param ret will be set to the intersection point (if any)
- * @return true if an intersection was found, false otherwise
+ * @return intersection if found, setNull() otherwise
*
*/
- bool lineIntersection(const LineSegment& line, Coordinate& ret) const;
+ Coordinate lineIntersection(const LineSegment& line) const;
/**
* Creates a LineString with the same coordinates as this segment
diff --git a/src/algorithm/CGAlgorithmsDD.cpp b/src/algorithm/CGAlgorithmsDD.cpp
index 7d884a9..c25c1f2 100644
--- a/src/algorithm/CGAlgorithmsDD.cpp
+++ b/src/algorithm/CGAlgorithmsDD.cpp
@@ -21,6 +21,7 @@
#include <geos/geom/Coordinate.h>
#include <geos/util/IllegalArgumentException.h>
#include <sstream>
+#include <cmath>
using namespace geos::geom;
using namespace geos::algorithm;
@@ -146,10 +147,9 @@ CGAlgorithmsDD::orientationIndexFilter(const Coordinate& pa,
return CGAlgorithmsDD::FAILURE;
}
-void
+Coordinate
CGAlgorithmsDD::intersection(const Coordinate& p1, const Coordinate& p2,
- const Coordinate& q1, const Coordinate& q2,
- Coordinate& rv)
+ const Coordinate& q1, const Coordinate& q2)
{
DD q1x(q1.x);
DD q1y(q1.y);
@@ -161,35 +161,32 @@ CGAlgorithmsDD::intersection(const Coordinate& p1, const Coordinate& p2,
DD p2x(p2.x);
DD p2y(p2.y);
- DD qdy = q2y - q1y;
- DD qdx = q2x - q1x;
+ DD px = p1y - p2y;
+ DD py = p2x - p1x;
+ DD pw = (p1x * p2y) - (p2x * p1y);
- DD pdy = p2y - p1y;
- DD pdx = p2x - p1x;
+ DD qx = q1y - q2y;
+ DD qy = q2x - q1x;
+ DD qw = (q1x * q2y) - (q2x * q1y);
- DD denom = (qdy * pdx) - (qdx * pdy);
+ DD x = (py * qw) - (qy * pw);
+ DD y = (qx * pw) - (px * qw);
+ DD w = (px * qy) - (qx * py);
- /**
- * Cases:
- * - denom is 0 if lines are parallel
- * - intersection point lies within line segment p if fracP is between 0 and 1
- * - intersection point lies within line segment q if fracQ is between 0 and 1
- */
- DD numx1 = (q2x - q1x) * (p1y - q1y);
- DD numx2 = (q2y - q1y) * (p1x - q1x);
- DD numx = numx1 - numx2;
- DD fracP = numx / denom;
+ double xInt = (x / w).ToDouble();
+ double yInt = (y / w).ToDouble();
- DD x = p1x + (p2x - p1x) * fracP;
+ Coordinate rv;
- DD numy1 = (p2x - p1x) * (p1y - q1y);
- DD numy2 = (p2y - p1y) * (p1x - q1x);
- DD numy = numy1 - numy2;
- DD fracQ = numy / denom;
- DD y = q1y + (q2y - q1y) * fracQ;
+ if (std::isnan(xInt) || std::isnan(yInt) ||
+ std::isinf(xInt) || std::isinf(yInt)) {
+ rv.setNull();
+ return rv;
+ }
- rv.x = x.ToDouble();
- rv.y = y.ToDouble();
+ rv.x = xInt;
+ rv.y = yInt;
+ return rv;
}
/* public static */
diff --git a/src/algorithm/Intersection.cpp b/src/algorithm/Intersection.cpp
new file mode 100644
index 0000000..834b7c0
--- /dev/null
+++ b/src/algorithm/Intersection.cpp
@@ -0,0 +1,88 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2019 Paul Ramsey <pramsey at cleverlephant.ca>
+ *
+ * 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.
+ *
+ **********************************************************************/
+
+#include <cmath>
+#include <vector>
+
+#include <geos/algorithm/Intersection.h>
+
+namespace geos {
+namespace algorithm { // geos.algorithm
+
+
+/* public static */
+geom::Coordinate
+Intersection::intersection(const geom::Coordinate& p1, const geom::Coordinate& p2,
+ const geom::Coordinate& q1, const geom::Coordinate& q2)
+{
+ double minX0 = p1.x < p2.x ? p1.x : p2.x;
+ double minY0 = p1.y < p2.y ? p1.y : p2.y;
+ double maxX0 = p1.x > p2.x ? p1.x : p2.x;
+ double maxY0 = p1.y > p2.y ? p1.y : p2.y;
+
+ double minX1 = q1.x < q2.x ? q1.x : q2.x;
+ double minY1 = q1.y < q2.y ? q1.y : q2.y;
+ double maxX1 = q1.x > q2.x ? q1.x : q2.x;
+ double maxY1 = q1.y > q2.y ? q1.y : q2.y;
+
+ double intMinX = minX0 > minX1 ? minX0 : minX1;
+ double intMaxX = maxX0 < maxX1 ? maxX0 : maxX1;
+ double intMinY = minY0 > minY1 ? minY0 : minY1;
+ double intMaxY = maxY0 < maxY1 ? maxY0 : maxY1;
+
+ double midx = (intMinX + intMaxX) / 2.0;
+ double midy = (intMinY + intMaxY) / 2.0;
+
+ // condition ordinate values by subtracting midpoint
+ double p1x = p1.x - midx;
+ double p1y = p1.y - midy;
+ double p2x = p2.x - midx;
+ double p2y = p2.y - midy;
+ double q1x = q1.x - midx;
+ double q1y = q1.y - midy;
+ double q2x = q2.x - midx;
+ double q2y = q2.y - midy;
+
+ // unrolled computation using homogeneous coordinates eqn
+ double px = p1y - p2y;
+ double py = p2x - p1x;
+ double pw = p1x * p2y - p2x * p1y;
+
+ double qx = q1y - q2y;
+ double qy = q2x - q1x;
+ double qw = q1x * q2y - q2x * q1y;
+
+ double x = py * qw - qy * pw;
+ double y = qx * pw - px * qw;
+ double w = px * qy - qx * py;
+
+ double xInt = x/w;
+ double yInt = y/w;
+ geom::Coordinate rv;
+ // check for parallel lines
+ if (std::isnan(xInt) || std::isnan(yInt) ||
+ std::isinf(xInt) || std::isinf(yInt)) {
+ rv.setNull();
+ return rv;
+ }
+ // de-condition intersection point
+ rv.x = xInt + midx;
+ rv.y = yInt + midy;
+ return rv;
+}
+
+
+} // namespace geos.algorithm
+} // namespace geos
+
diff --git a/src/algorithm/LineIntersector.cpp b/src/algorithm/LineIntersector.cpp
index e928e5b..33f379d 100644
--- a/src/algorithm/LineIntersector.cpp
+++ b/src/algorithm/LineIntersector.cpp
@@ -21,7 +21,7 @@
#include <geos/algorithm/LineIntersector.h>
#include <geos/algorithm/Distance.h>
#include <geos/algorithm/Orientation.h>
-#include <geos/algorithm/HCoordinate.h>
+#include <geos/algorithm/Intersection.h>
#include <geos/algorithm/NotRepresentableException.h>
#include <geos/geom/Coordinate.h>
#include <geos/geom/PrecisionModel.h>
@@ -385,8 +385,8 @@ LineIntersector::hasIntersection(const Coordinate& p, const Coordinate& p1, cons
/*private*/
int
-LineIntersector::computeIntersect(const Coordinate& p1, const Coordinate& p2, const Coordinate& q1,
- const Coordinate& q2)
+LineIntersector::computeIntersect(const Coordinate& p1, const Coordinate& p2,
+ const Coordinate& q1, const Coordinate& q2)
{
#if GEOS_DEBUG
cerr << "LineIntersector::computeIntersect called" << endl;
@@ -543,7 +543,7 @@ LineIntersector::computeIntersect(const Coordinate& p1, const Coordinate& p2, co
}
else {
isProperVar = true;
- intersection(p1, p2, q1, q2, intPt[0]);
+ intPt[0] = intersection(p1, p2, q1, q2);
}
#if GEOS_DEBUG
cerr << " POINT_INTERSECTION; intPt[0]:" << intPt[0].toString() << endl;
@@ -553,8 +553,8 @@ LineIntersector::computeIntersect(const Coordinate& p1, const Coordinate& p2, co
/*private*/
int
-LineIntersector::computeCollinearIntersection(const Coordinate& p1, const Coordinate& p2, const Coordinate& q1,
- const Coordinate& q2)
+LineIntersector::computeCollinearIntersection(const Coordinate& p1, const Coordinate& p2,
+ const Coordinate& q1, const Coordinate& q2)
{
#if COMPUTE_Z
double ztot;
@@ -840,13 +840,12 @@ LineIntersector::computeCollinearIntersection(const Coordinate& p1, const Coordi
}
/*private*/
-void
-LineIntersector::intersection(const Coordinate& p1,
- const Coordinate& p2, const Coordinate& q1, const Coordinate& q2,
- Coordinate& intPtOut) const
+Coordinate
+LineIntersector::intersection(const Coordinate& p1, const Coordinate& p2,
+ const Coordinate& q1, const Coordinate& q2) const
{
- intersectionWithNormalization(p1, p2, q1, q2, intPtOut);
+ Coordinate intPtOut = intersectionSafe(p1, p2, q1, q2);
/*
* Due to rounding it can happen that the computed intersection is
@@ -894,47 +893,7 @@ LineIntersector::intersection(const Coordinate& p1,
intPtOut.z = ztot / zvals;
}
#endif // COMPUTE_Z
-
-}
-
-/*private*/
-void
-LineIntersector::intersectionWithNormalization(const Coordinate& p1,
- const Coordinate& p2, const Coordinate& q1, const Coordinate& q2,
- Coordinate& intPtOut) const
-{
- Coordinate n1 = p1;
- Coordinate n2 = p2;
- Coordinate n3 = q1;
- Coordinate n4 = q2;
- Coordinate normPt;
- normalizeToEnvCentre(n1, n2, n3, n4, normPt);
-
- safeHCoordinateIntersection(n1, n2, n3, n4, intPtOut);
-
- intPtOut.x += normPt.x;
- intPtOut.y += normPt.y;
-}
-
-
-/*private*/
-double
-LineIntersector::smallestInAbsValue(double x1, double x2, double x3, double x4) const
-{
- double x = x1;
- double xabs = fabs(x);
- if(fabs(x2) < xabs) {
- x = x2;
- xabs = fabs(x2);
- }
- if(fabs(x3) < xabs) {
- x = x3;
- xabs = fabs(x3);
- }
- if(fabs(x4) < xabs) {
- x = x4;
- }
- return x;
+ return intPtOut;
}
/*private*/
@@ -947,82 +906,17 @@ LineIntersector::isInSegmentEnvelopes(const Coordinate& pt) const
}
/*private*/
-void
-LineIntersector::normalizeToEnvCentre(Coordinate& n00, Coordinate& n01,
- Coordinate& n10, Coordinate& n11, Coordinate& normPt) const
+Coordinate
+LineIntersector::intersectionSafe(const Coordinate& p1, const Coordinate& p2,
+ const Coordinate& q1, const Coordinate& q2) const
{
- double minX0 = n00.x < n01.x ? n00.x : n01.x;
- double minY0 = n00.y < n01.y ? n00.y : n01.y;
- double maxX0 = n00.x > n01.x ? n00.x : n01.x;
- double maxY0 = n00.y > n01.y ? n00.y : n01.y;
-
- double minX1 = n10.x < n11.x ? n10.x : n11.x;
- double minY1 = n10.y < n11.y ? n10.y : n11.y;
- double maxX1 = n10.x > n11.x ? n10.x : n11.x;
- double maxY1 = n10.y > n11.y ? n10.y : n11.y;
-
- double intMinX = minX0 > minX1 ? minX0 : minX1;
- double intMaxX = maxX0 < maxX1 ? maxX0 : maxX1;
- double intMinY = minY0 > minY1 ? minY0 : minY1;
- double intMaxY = maxY0 < maxY1 ? maxY0 : maxY1;
-
- double intMidX = (intMinX + intMaxX) / 2.0;
- double intMidY = (intMinY + intMaxY) / 2.0;
-
- normPt.x = intMidX;
- normPt.y = intMidY;
-
- n00.x -= normPt.x;
- n00.y -= normPt.y;
- n01.x -= normPt.x;
- n01.y -= normPt.y;
- n10.x -= normPt.x;
- n10.y -= normPt.y;
- n11.x -= normPt.x;
- n11.y -= normPt.y;
-
-#if COMPUTE_Z
-
- // Only do this if input does have Z
- // See https://trac.osgeo.org/geos/ticket/811
- if(std::isnan(n00.z)) {
- return;
+ Coordinate intPt = Intersection::intersection(p1, p2, q1, q2);
+ if (intPt.isNull()) {
+ intPt = nearestEndpoint(p1, p2, q1, q2);
}
-
- double minZ0 = n00.z < n01.z ? n00.z : n01.z;
- double minZ1 = n10.z < n11.z ? n10.z : n11.z;
- double maxZ0 = n00.z > n01.z ? n00.z : n01.z;
- double maxZ1 = n10.z > n11.z ? n10.z : n11.z;
- double intMinZ = minZ0 > minZ1 ? minZ0 : minZ1;
- double intMaxZ = maxZ0 < maxZ1 ? maxZ0 : maxZ1;
- double intMidZ = (intMinZ + intMaxZ) / 2.0;
- normPt.z = intMidZ;
- n00.z -= normPt.z;
- n01.z -= normPt.z;
- n10.z -= normPt.z;
- n11.z -= normPt.z;
-#endif
+ return intPt;
}
-/*private*/
-void
-LineIntersector::safeHCoordinateIntersection(const Coordinate& p1,
- const Coordinate& p2, const Coordinate& q1,
- const Coordinate& q2, Coordinate& intPtOut) const
-{
- try {
- HCoordinate::intersection(p1, p2, q1, q2, intPtOut);
-#if GEOS_DEBUG
- cerr << " HCoordinate found intersection h:" << intPt.toString() << endl;
-#endif
-
- }
- catch(const NotRepresentableException& /* e */) {
- // compute an approximate result
- //intPt = CentralEndpointIntersector::getIntersection(p1, p2, q1, q2);
- intPtOut = nearestEndpoint(p1, p2, q1, q2);
- }
-}
} // namespace geos.algorithm
} // namespace geos
diff --git a/src/algorithm/Makefile.am b/src/algorithm/Makefile.am
index 2579281..75dfb6c 100644
--- a/src/algorithm/Makefile.am
+++ b/src/algorithm/Makefile.am
@@ -21,6 +21,7 @@ libalgorithm_la_SOURCES = \
InteriorPointArea.cpp \
InteriorPointLine.cpp \
InteriorPointPoint.cpp \
+ Intersection.cpp \
LineIntersector.cpp \
Length.cpp \
MinimumBoundingCircle.cpp \
@@ -33,7 +34,6 @@ libalgorithm_la_SOURCES = \
RayCrossingCounterDD.cpp \
RobustDeterminant.cpp
-
libalgorithm_la_LIBADD = \
locate/liblocation.la \
distance/libdistance.la
diff --git a/src/algorithm/MinimumDiameter.cpp b/src/algorithm/MinimumDiameter.cpp
index 762d3a9..31ceb24 100644
--- a/src/algorithm/MinimumDiameter.cpp
+++ b/src/algorithm/MinimumDiameter.cpp
@@ -343,11 +343,10 @@ MinimumDiameter::getMinimumRectangle()
LineSegment minParaLine = computeSegmentForLine(-dy, dx, minPara);
// compute vertices of rectangle (where the para/perp max & min lines intersect)
- Coordinate p0, p1, p2, p3;
- maxParaLine.lineIntersection(maxPerpLine, p0);
- minParaLine.lineIntersection(maxPerpLine, p1);
- minParaLine.lineIntersection(minPerpLine, p2);
- maxParaLine.lineIntersection(minPerpLine, p3);
+ Coordinate p0 = maxParaLine.lineIntersection(maxPerpLine);
+ Coordinate p1 = minParaLine.lineIntersection(maxPerpLine);
+ Coordinate p2 = minParaLine.lineIntersection(minPerpLine);
+ Coordinate p3 = maxParaLine.lineIntersection(minPerpLine);
const CoordinateSequenceFactory* csf =
inputGeom->getFactory()->getCoordinateSequenceFactory();
diff --git a/src/geom/LineSegment.cpp b/src/geom/LineSegment.cpp
index 1b6bb5d..699316f 100644
--- a/src/geom/LineSegment.cpp
+++ b/src/geom/LineSegment.cpp
@@ -27,8 +27,7 @@
#include <geos/geom/CoordinateArraySequence.h> // should we really be using this?
#include <geos/algorithm/Orientation.h>
#include <geos/algorithm/LineIntersector.h>
-#include <geos/algorithm/HCoordinate.h>
-#include <geos/algorithm/NotRepresentableException.h>
+#include <geos/algorithm/Intersection.h>
#include <geos/util/IllegalStateException.h>
#include <geos/profiler.h>
#include <geos/inline.h>
@@ -41,11 +40,6 @@
# include <geos/geom/LineSegment.inl>
#endif
-using namespace std;
-//using namespace geos::algorithm;
-using geos::algorithm::HCoordinate;
-using geos::algorithm::NotRepresentableException;
-using geos::algorithm::LineIntersector;
namespace geos {
namespace geom { // geos::geom
@@ -194,8 +188,8 @@ std::array<Coordinate, 2>
LineSegment::closestPoints(const LineSegment& line)
{
// test for intersection
- Coordinate intPt;
- if(intersection(line, intPt)) {
+ Coordinate intPt = intersection(line);
+ if(!intPt.isNull()) {
return { intPt, intPt };
}
@@ -245,29 +239,23 @@ LineSegment::closestPoints(const LineSegment& line)
return closestPt;
}
-bool
-LineSegment::intersection(const LineSegment& line, Coordinate& ret) const
+Coordinate
+LineSegment::intersection(const LineSegment& line) const
{
algorithm::LineIntersector li;
li.computeIntersection(p0, p1, line.p0, line.p1);
if(li.hasIntersection()) {
- ret = li.getIntersection(0);
- return true;
+ return li.getIntersection(0);
}
- return false;
+ Coordinate rv;
+ rv.setNull();
+ return rv;
}
-bool
-LineSegment::lineIntersection(const LineSegment& line, Coordinate& ret) const
+Coordinate
+LineSegment::lineIntersection(const LineSegment& line) const
{
- try {
- HCoordinate::intersection(p0, p1, line.p0, line.p1, ret);
- return true;
- }
- catch(const NotRepresentableException& /*ex*/) {
- // eat this exception, and return null;
- }
- return false;
+ return algorithm::Intersection::intersection(p0, p1, line.p0, line.p1);
}
diff --git a/src/operation/buffer/OffsetSegmentGenerator.cpp b/src/operation/buffer/OffsetSegmentGenerator.cpp
index 6c392ae..9605ec1 100644
--- a/src/operation/buffer/OffsetSegmentGenerator.cpp
+++ b/src/operation/buffer/OffsetSegmentGenerator.cpp
@@ -33,7 +33,7 @@
#include <geos/geom/Coordinate.h>
#include <geos/geom/PrecisionModel.h>
#include <geos/algorithm/NotRepresentableException.h>
-#include <geos/algorithm/HCoordinate.h>
+#include <geos/algorithm/Intersection.h>
#include <geos/util.h>
#ifndef GEOS_DEBUG
@@ -481,43 +481,22 @@ OffsetSegmentGenerator::addMitreJoin(const geom::Coordinate& p,
const geom::LineSegment& p_offset1,
double p_distance)
{
- bool isMitreWithinLimit = true;
- Coordinate intPt;
-
/**
- * This computation is unstable if the offset segments
- * are nearly collinear.
- * Howver, this situation should have been eliminated earlier
- * by the check for whether the offset segment endpoints are
- * almost coincident
+ * This computation is unstable if the offset segments are nearly collinear.
+ * However, this situation should have been eliminated earlier by the check
+ * for whether the offset segment endpoints are almost coincident
*/
- try {
- HCoordinate::intersection(p_offset0.p0, p_offset0.p1,
- p_offset1.p0, p_offset1.p1,
- intPt);
-
- double mitreRatio = p_distance <= 0.0 ? 1.0
- : intPt.distance(p) / fabs(p_distance);
-
- if(mitreRatio > bufParams.getMitreLimit()) {
- isMitreWithinLimit = false;
+ Coordinate intPt = algorithm::Intersection::intersection(offset0.p0, offset0.p1, offset1.p0, offset1.p1);
+ if (!intPt.isNull()) {
+ double mitreRatio = p_distance <= 0.0 ? 1.0 : intPt.distance(p) / fabs(p_distance);
+ if (mitreRatio <= bufParams.getMitreLimit()) {
+ segList.addPt(intPt);
+ return;
}
}
- catch(const NotRepresentableException& e) {
- ::geos::ignore_unused_variable_warning(e);
-
- intPt = Coordinate(0, 0);
- isMitreWithinLimit = false;
- }
-
- if(isMitreWithinLimit) {
- segList.addPt(intPt);
- }
- else {
- addLimitedMitreJoin(p_offset0, p_offset1, p_distance,
- bufParams.getMitreLimit());
- //addBevelJoin(offset0, offset1);
- }
+ // at this point either intersection failed or mitre limit was exceeded
+ addLimitedMitreJoin(p_offset0, p_offset1, p_distance,
+ bufParams.getMitreLimit());
}
/* private */
diff --git a/tests/unit/Makefile.am b/tests/unit/Makefile.am
index 71e1210..a7184de 100644
--- a/tests/unit/Makefile.am
+++ b/tests/unit/Makefile.am
@@ -49,6 +49,7 @@ geos_unit_SOURCES = \
algorithm/distance/DiscreteFrechetDistanceTest.cpp \
algorithm/distance/DiscreteHausdorffDistanceTest.cpp \
algorithm/InteriorPointAreaTest.cpp \
+ algorithm/IntersectionTest.cpp \
algorithm/LengthTest.cpp \
algorithm/LocatePointInRingTest.cpp \
algorithm/MinimumBoundingCircleTest.cpp \
diff --git a/tests/unit/algorithm/IntersectionTest.cpp b/tests/unit/algorithm/IntersectionTest.cpp
new file mode 100644
index 0000000..9cae8e4
--- /dev/null
+++ b/tests/unit/algorithm/IntersectionTest.cpp
@@ -0,0 +1,117 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2019 Martin Davis <mtnclimb at gmail.com>
+ * Copyright (C) 2011 Sandro Santilli <strk at kbt.io>
+ *
+ * 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.
+ *
+ **********************************************************************/
+
+//
+// Test Suite for geos::algorithm::InteriorPointArea
+
+#include <tut/tut.hpp>
+// geos
+#include <geos/geom/Coordinate.h>
+#include <geos/algorithm/Intersection.h>
+#include <geos/io/WKTReader.h>
+#include <geos/geom/Geometry.h>
+// std
+#include <sstream>
+#include <string>
+#include <memory>
+
+namespace tut {
+//
+// Test Group
+//
+
+// dummy data, not used
+struct test_intersection_data {
+ typedef geos::geom::Geometry Geometry;
+ typedef geos::geom::Coordinate Coordinate;
+ typedef geos::algorithm::Intersection Intersection;
+
+ geos::io::WKTReader reader;
+ std::unique_ptr<Geometry> geom;
+
+ double MAX_ABS_ERROR = 1e-5;
+
+ 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);
+ Coordinate q2(q2x, q2y);
+ Coordinate actual = Intersection::intersection(p1, p2, q1, q2);
+ 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) {
+ Coordinate p1(p1x, p1y);
+ Coordinate p2(p2x, p2y);
+ Coordinate q1(q1x, q1y);
+ Coordinate q2(q2x, q2y);
+ Coordinate expected(expectedx, expectedy);
+ Coordinate actual = Intersection::intersection(p1, p2, q1, q2);
+ double dist = actual.distance(expected);
+ // std::cout << "Expected: " << expected << " Actual: " << actual << " Dist = " << dist << std::endl;
+ ensure("checkIntersection", dist <= MAX_ABS_ERROR);
+ }
+
+ test_intersection_data()
+ {}
+
+};
+
+typedef test_group<test_intersection_data> group;
+typedef group::object object;
+
+group test_intersection_data("geos::algorithm::Intersection");
+
+
+template<>
+template<>
+void object::test<1>
+()
+{
+ // testSimple
+ checkIntersection(
+ 0,0, 10,10,
+ 0,10, 10,0,
+ 5,5);
+
+ // testCollinear
+ checkIntersectionNull(
+ 0,0, 10,10,
+ 20,20, 30, 30 );
+
+ // testParallel
+ checkIntersectionNull(
+ 0,0, 10,10,
+ 10,0, 20,10 );
+
+ // testAlmostCollinear
+ checkIntersection(
+ 35613471.6165017, 4257145.306132293, 35613477.7705378, 4257160.528222711,
+ 35613477.77505724, 4257160.539653536, 35613479.85607389, 4257165.92369170,
+ 35613477.772841461, 4257160.5339209242 );
+
+ // testAlmostCollinearCond
+ checkIntersection(
+ 1.6165017, 45.306132293, 7.7705378, 60.528222711,
+ 7.77505724, 60.539653536, 9.85607389, 65.92369170,
+ 7.772841461, 60.5339209242 );
+
+}
+
+} // namespace tut
+
diff --git a/tests/unit/geom/LineSegmentTest.cpp b/tests/unit/geom/LineSegmentTest.cpp
index 701d2e2..eb95aae 100644
--- a/tests/unit/geom/LineSegmentTest.cpp
+++ b/tests/unit/geom/LineSegmentTest.cpp
@@ -14,12 +14,30 @@ namespace tut {
//
struct test_lineseg_data {
+
+ typedef geos::geom::Coordinate Coordinate;
+ typedef geos::geom::LineSegment LineSegment;
+
geos::geom::Coordinate ph1;
geos::geom::Coordinate ph2;
geos::geom::Coordinate pv1;
geos::geom::Coordinate pv2;
geos::geom::LineSegment h1;
geos::geom::LineSegment v1;
+ double MAX_ABS_ERROR_INTERSECTION = 1e-5;
+
+ void checkLineIntersection(double p1x, double p1y, double p2x, double p2y,
+ double q1x, double q1y, double q2x, double q2y,
+ double expectedx, double expectedy) {
+ LineSegment seg1(p1x, p1y, p2x, p2y);
+ LineSegment seg2(q1x, q1y, q2x, q2y);
+
+ Coordinate actual = seg1.lineIntersection(seg2);
+ Coordinate expected(expectedx, expectedy);
+ double dist = actual.distance(expected);
+ // std::cout << "Expected: " << expected << " Actual: " << actual << " Dist = " << dist << std::endl;
+ ensure("checkLineIntersection", dist <= MAX_ABS_ERROR_INTERSECTION);
+ }
test_lineseg_data()
: ph1(0, 2), ph2(10, 2), pv1(0, 0), pv2(0, 10), h1(ph1, ph2), v1(pv1, pv2)
@@ -117,5 +135,23 @@ void object::test<6>
ensure_equals(v1.distance(p), 1);
}
+template<>
+template<>
+void object::test<7>
+()
+{
+ // simple case
+ checkLineIntersection(
+ 0,0, 10,10,
+ 0,10, 10,0,
+ 5,5);
+
+ //Almost collinear - See JTS GitHub issue #464
+ checkLineIntersection(
+ 35613471.6165017, 4257145.306132293, 35613477.7705378, 4257160.528222711,
+ 35613477.77505724, 4257160.539653536, 35613479.85607389, 4257165.92369170,
+ 35613477.772841461, 4257160.5339209242 );
+}
+
} // namespace tut
-----------------------------------------------------------------------
Summary of changes:
capi/geos_ts_c.cpp | 6 +-
include/geos/algorithm/CGAlgorithmsDD.h | 14 ++-
include/geos/algorithm/Intersection.h | 61 ++++++++++
include/geos/algorithm/LineIntersector.h | 52 ++-------
include/geos/algorithm/Makefile.am | 1 +
include/geos/geom/LineSegment.h | 10 +-
src/algorithm/CGAlgorithmsDD.cpp | 49 ++++----
src/algorithm/Intersection.cpp | 88 +++++++++++++++
src/algorithm/LineIntersector.cpp | 142 +++---------------------
src/algorithm/Makefile.am | 2 +-
src/algorithm/MinimumDiameter.cpp | 9 +-
src/geom/LineSegment.cpp | 36 ++----
src/operation/buffer/OffsetSegmentGenerator.cpp | 47 +++-----
tests/unit/Makefile.am | 1 +
tests/unit/algorithm/IntersectionTest.cpp | 117 +++++++++++++++++++
tests/unit/geom/LineSegmentTest.cpp | 36 ++++++
16 files changed, 404 insertions(+), 267 deletions(-)
create mode 100644 include/geos/algorithm/Intersection.h
create mode 100644 src/algorithm/Intersection.cpp
create mode 100644 tests/unit/algorithm/IntersectionTest.cpp
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list