[geos-commits] [SCM] GEOS branch master updated. 03e8090d9a2019402dbb0d3c2931eac75d11030e
git at osgeo.org
git at osgeo.org
Mon Nov 9 01:14:18 PST 2020
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 03e8090d9a2019402dbb0d3c2931eac75d11030e (commit)
from 26cd4782faf7451c7a3df6a33a462b2669ef5153 (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 03e8090d9a2019402dbb0d3c2931eac75d11030e
Author: Sandro Santilli <strk at kbt.io>
Date: Sat Nov 7 12:19:44 2020 +0100
Port Z handling in JTS Intersection
References #1063
diff --git a/include/geos/algorithm/LineIntersector.h b/include/geos/algorithm/LineIntersector.h
index 8b1232e..2ac07fe 100644
--- a/include/geos/algorithm/LineIntersector.h
+++ b/include/geos/algorithm/LineIntersector.h
@@ -343,8 +343,56 @@ private:
geom::Coordinate intersectionSafe(const geom::Coordinate& p1, const geom::Coordinate& p2,
const geom::Coordinate& q1, const geom::Coordinate& q2) const;
+ /**
+ * Finds the endpoint of the segments P and Q which
+ * is closest to the other segment.
+ * This is a reasonable surrogate for the true
+ * intersection points in ill-conditioned cases
+ * (e.g. where two segments are nearly coincident,
+ * or where the endpoint of one segment lies almost on the other segment).
+ * <p>
+ * This replaces the older CentralEndpoint heuristic,
+ * which chose the wrong endpoint in some cases
+ * where the segments had very distinct slopes
+ * and one endpoint lay almost on the other segment.
+ *
+ * @param p1 an endpoint of segment P
+ * @param p2 an endpoint of segment P
+ * @param q1 an endpoint of segment Q
+ * @param q2 an endpoint of segment Q
+ * @return the nearest endpoint to the other segment
+ */
+ static geom::Coordinate nearestEndpoint(const geom::Coordinate& p1,
+ const geom::Coordinate& p2,
+ const geom::Coordinate& q1,
+ const geom::Coordinate& q2);
+
+ static double zGet(const geom::Coordinate& p, const geom::Coordinate& q);
+
+ static double zGetOrInterpolate(const geom::Coordinate& p,
+ const geom::Coordinate& p0,
+ const geom::Coordinate& p1);
+
+ static geom::Coordinate zGetOrInterpolateCopy(const geom::Coordinate& p,
+ const geom::Coordinate& p0,
+ const geom::Coordinate& p1);
+
+ /// \brief
+ /// Return a Z value being the interpolation of Z from p0 to p1 at
+ /// the given point p
+ static double zInterpolate(const geom::Coordinate& p,
+ const geom::Coordinate& p0,
+ const geom::Coordinate& p1);
+
+ static double zInterpolate(const geom::Coordinate& p,
+ const geom::Coordinate& p1,
+ const geom::Coordinate& p2,
+ const geom::Coordinate& q1,
+ const geom::Coordinate& q2);
+
};
+
} // namespace geos::algorithm
} // namespace geos
diff --git a/src/algorithm/LineIntersector.cpp b/src/algorithm/LineIntersector.cpp
index 10878b7..080590c 100644
--- a/src/algorithm/LineIntersector.cpp
+++ b/src/algorithm/LineIntersector.cpp
@@ -41,10 +41,6 @@
#include <iostream>
#endif
-#ifndef COMPUTE_Z
-#define COMPUTE_Z 1
-#endif // COMPUTE_Z
-
using namespace std;
using namespace geos::geom;
@@ -52,54 +48,6 @@ using namespace geos::geom;
namespace geos {
namespace algorithm { // geos.algorithm
-namespace { // anonymous
-
-/**
- * Finds the endpoint of the segments P and Q which
- * is closest to the other segment.
- * This is a reasonable surrogate for the true
- * intersection points in ill-conditioned cases
- * (e.g. where two segments are nearly coincident,
- * or where the endpoint of one segment lies almost on the other segment).
- * <p>
- * This replaces the older CentralEndpoint heuristic,
- * which chose the wrong endpoint in some cases
- * where the segments had very distinct slopes
- * and one endpoint lay almost on the other segment.
- *
- * @param p1 an endpoint of segment P
- * @param p2 an endpoint of segment P
- * @param q1 an endpoint of segment Q
- * @param q2 an endpoint of segment Q
- * @return the nearest endpoint to the other segment
- */
-Coordinate
-nearestEndpoint(const Coordinate& p1, const Coordinate& p2,
- const Coordinate& q1, const Coordinate& q2)
-{
- const Coordinate* nearestPt = &p1;
- double minDist = Distance::pointToSegment(p1, q1, q2);
-
- double dist = Distance::pointToSegment(p2, q1, q2);
- if(dist < minDist) {
- minDist = dist;
- nearestPt = &p2;
- }
- dist = Distance::pointToSegment(q1, p1, p2);
- if(dist < minDist) {
- minDist = dist;
- nearestPt = &q1;
- }
- dist = Distance::pointToSegment(q2, p1, p2);
- if(dist < minDist) {
- nearestPt = &q2;
- }
- return *nearestPt;
-}
-
-
-} // anonymous namespace
-
/*public static*/
double
LineIntersector::computeEdgeDistance(const Coordinate& p, const Coordinate& p0, const Coordinate& p1)
@@ -347,21 +295,6 @@ LineIntersector::computeIntersection(const Coordinate& p, const Coordinate& p1,
if((p == p1) || (p == p2)) { // 2d only test
isProperVar = false;
}
-#if COMPUTE_Z
- intPt[0] = p;
-#if GEOS_DEBUG
- cerr << "RobustIntersector::computeIntersection(Coordinate,Coordinate,Coordinate) calling interpolateZ" << endl;
-#endif
- double z = interpolateZ(p, p1, p2);
- if(!std::isnan(z)) {
- if(std::isnan(intPt[0].z)) {
- intPt[0].z = z;
- }
- else {
- intPt[0].z = (intPt[0].z + z) / 2;
- }
- }
-#endif // COMPUTE_Z
result = POINT_INTERSECTION;
return;
}
@@ -426,6 +359,9 @@ LineIntersector::computeIntersect(const Coordinate& p1, const Coordinate& p2,
return NO_INTERSECTION;
}
+ /**
+ * Intersection is collinear if each endpoint lies on the other line.
+ */
bool collinear = Pq1 == 0 && Pq2 == 0 && Qp1 == 0 && Qp2 == 0;
if(collinear) {
#if GEOS_DEBUG
@@ -449,13 +385,17 @@ LineIntersector::computeIntersect(const Coordinate& p1, const Coordinate& p2,
* since at this point we know that the inputLines must
* intersect.
*/
+ Coordinate p;
+ double z = std::numeric_limits<double>::quiet_NaN();
+
if(Pq1 == 0 || Pq2 == 0 || Qp1 == 0 || Qp2 == 0) {
-#if COMPUTE_Z
- int hits = 0;
- double z = 0.0;
-#endif
+
isProperVar = false;
+#if GEOS_DEBUG
+ cerr << " intersection is NOT proper" << endl;
+#endif
+
/* Check for two equal endpoints.
* This is done explicitly rather than by the orientation tests
* below in order to improve robustness.
@@ -469,81 +409,61 @@ LineIntersector::computeIntersect(const Coordinate& p1, const Coordinate& p2,
* LINESTRING ( -48.51001596420236 -22.063180333403878,
* 19.850257749638203 46.29709338043669 )
*
- * which used to produce the result:
+ * which used to produce the INCORRECT result:
* (20.31970698357233, 46.76654261437082, NaN)
*/
- if(p1.equals2D(q1) || p1.equals2D(q2)) {
- intPt[0] = p1;
-#if COMPUTE_Z
- if(!std::isnan(p1.z)) {
- z += p1.z;
- hits++;
- }
-#endif
+ if (p1.equals2D(q1)) {
+ p = p1;
+ z = zGet(p1, q1);
}
- else if(p2.equals2D(q1) || p2.equals2D(q2)) {
- intPt[0] = p2;
-#if COMPUTE_Z
- if(!std::isnan(p2.z)) {
- z += p2.z;
- hits++;
- }
-#endif
+ else if (p1.equals2D(q2)) {
+ p = p1;
+ z = zGet(p1, q2);
+ }
+ else if (p2.equals2D(q1)) {
+ p = p2;
+ z = zGet(p2, q1);
+ }
+ else if (p2.equals2D(q2)) {
+ p = p2;
+ z = zGet(p2, q2);
}
-
/*
* Now check to see if any endpoint lies on the interior of the other segment.
*/
else if(Pq1 == 0) {
- intPt[0] = q1;
-#if COMPUTE_Z
- if(!std::isnan(q1.z)) {
- z += q1.z;
- hits++;
- }
-#endif
+ p = q1;
+ z = zGetOrInterpolate(q1, p1, p2);
}
else if(Pq2 == 0) {
- intPt[0] = q2;
-#if COMPUTE_Z
- if(!std::isnan(q2.z)) {
- z += q2.z;
- hits++;
- }
-#endif
+ p = q2;
+ z = zGetOrInterpolate(q2, p1, p2);
}
else if(Qp1 == 0) {
- intPt[0] = p1;
-#if COMPUTE_Z
- if(!std::isnan(p1.z)) {
- z += p1.z;
- hits++;
- }
-#endif
+ p = p1;
+ z = zGetOrInterpolate(p1, q1, q2);
}
else if(Qp2 == 0) {
- intPt[0] = p2;
-#if COMPUTE_Z
- if(!std::isnan(p2.z)) {
- z += p2.z;
- hits++;
- }
-#endif
+ p = p2;
+ z = zGetOrInterpolate(p2, q1, q2);
}
-#if COMPUTE_Z
-#if GEOS_DEBUG
- cerr << "LineIntersector::computeIntersect: z:" << z << " hits:" << hits << endl;
-#endif // GEOS_DEBUG
- if(hits) {
- intPt[0].z = z / hits;
- }
-#endif // COMPUTE_Z
}
else {
+#if GEOS_DEBUG
+ cerr << " intersection is proper" << endl;
+#endif // GEOS_DEBUG
isProperVar = true;
- intPt[0] = intersection(p1, p2, q1, q2);
+ p = intersection(p1, p2, q1, q2);
+#if GEOS_DEBUG
+ cerr << " computed intersection point: " << p << endl;
+#endif // GEOS_DEBUG
+ z = zInterpolate(p, p1, p2, q1, q2);
+#if GEOS_DEBUG
+ cerr << " computed proper intersection Z: " << z << endl;
+#endif // GEOS_DEBUG
}
+ intPt[0] = Coordinate(p.x, p.y, z);
#if GEOS_DEBUG
cerr << " POINT_INTERSECTION; intPt[0]:" << intPt[0].toString() << endl;
#endif // GEOS_DEBUG
@@ -555,14 +475,6 @@ int
LineIntersector::computeCollinearIntersection(const Coordinate& p1, const Coordinate& p2,
const Coordinate& q1, const Coordinate& q2)
{
-#if COMPUTE_Z
- double ztot;
- int hits;
- double p2z;
- double p1z;
- double q1z;
- double q2z;
-#endif // COMPUTE_Z
#if GEOS_DEBUG
cerr << "LineIntersector::computeCollinearIntersection called" << endl;
@@ -570,270 +482,84 @@ LineIntersector::computeCollinearIntersection(const Coordinate& p1, const Coordi
endl;
#endif // GEOS_DEBUG
- bool p1q1p2 = Envelope::intersects(p1, p2, q1);
- bool p1q2p2 = Envelope::intersects(p1, p2, q2);
- bool q1p1q2 = Envelope::intersects(q1, q2, p1);
- bool q1p2q2 = Envelope::intersects(q1, q2, p2);
+ bool q1inP = Envelope::intersects(p1, p2, q1);
+ bool q2inP = Envelope::intersects(p1, p2, q2);
+ bool p1inQ = Envelope::intersects(q1, q2, p1);
+ bool p2inQ = Envelope::intersects(q1, q2, p2);
- if(p1q1p2 && p1q2p2) {
+ if(q1inP && q2inP) {
#if GEOS_DEBUG
- cerr << " p1q1p2 && p1q2p2" << endl;
-#endif
- intPt[0] = q1;
-#if COMPUTE_Z
- ztot = 0;
- hits = 0;
- q1z = interpolateZ(q1, p1, p2);
- if(!std::isnan(q1z)) {
- ztot += q1z;
- hits++;
- }
- if(!std::isnan(q1.z)) {
- ztot += q1.z;
- hits++;
- }
- if(hits) {
- intPt[0].z = ztot / hits;
- }
-#endif
- intPt[1] = q2;
-#if COMPUTE_Z
- ztot = 0;
- hits = 0;
- q2z = interpolateZ(q2, p1, p2);
- if(!std::isnan(q2z)) {
- ztot += q2z;
- hits++;
- }
- if(!std::isnan(q2.z)) {
- ztot += q2.z;
- hits++;
- }
- if(hits) {
- intPt[1].z = ztot / hits;
- }
+ cerr << " q1inP && q2inP" << endl;
#endif
+ intPt[0] = zGetOrInterpolateCopy(q1, p1, p2);
+ intPt[1] = zGetOrInterpolateCopy(q2, p1, p2);
#if GEOS_DEBUG
cerr << " intPt[0]: " << intPt[0].toString() << endl;
cerr << " intPt[1]: " << intPt[1].toString() << endl;
#endif
return COLLINEAR_INTERSECTION;
}
- if(q1p1q2 && q1p2q2) {
+ if(p1inQ && p2inQ) {
#if GEOS_DEBUG
- cerr << " q1p1q2 && q1p2q2" << endl;
-#endif
- intPt[0] = p1;
-#if COMPUTE_Z
- ztot = 0;
- hits = 0;
- p1z = interpolateZ(p1, q1, q2);
- if(!std::isnan(p1z)) {
- ztot += p1z;
- hits++;
- }
- if(!std::isnan(p1.z)) {
- ztot += p1.z;
- hits++;
- }
- if(hits) {
- intPt[0].z = ztot / hits;
- }
-#endif
- intPt[1] = p2;
-#if COMPUTE_Z
- ztot = 0;
- hits = 0;
- p2z = interpolateZ(p2, q1, q2);
- if(!std::isnan(p2z)) {
- ztot += p2z;
- hits++;
- }
- if(!std::isnan(p2.z)) {
- ztot += p2.z;
- hits++;
- }
- if(hits) {
- intPt[1].z = ztot / hits;
- }
+ cerr << " p1inQ && p2inQ" << endl;
#endif
+ intPt[0] = zGetOrInterpolateCopy(p1, q1, q2);
+ intPt[1] = zGetOrInterpolateCopy(p2, q1, q2);
return COLLINEAR_INTERSECTION;
}
- if(p1q1p2 && q1p1q2) {
+ if(q1inP && p1inQ) {
#if GEOS_DEBUG
- cerr << " p1q1p2 && q1p1q2" << endl;
-#endif
- intPt[0] = q1;
-#if COMPUTE_Z
- ztot = 0;
- hits = 0;
- q1z = interpolateZ(q1, p1, p2);
- if(!std::isnan(q1z)) {
- ztot += q1z;
- hits++;
- }
- if(!std::isnan(q1.z)) {
- ztot += q1.z;
- hits++;
- }
- if(hits) {
- intPt[0].z = ztot / hits;
- }
-#endif
- intPt[1] = p1;
-#if COMPUTE_Z
- ztot = 0;
- hits = 0;
- p1z = interpolateZ(p1, q1, q2);
- if(!std::isnan(p1z)) {
- ztot += p1z;
- hits++;
- }
- if(!std::isnan(p1.z)) {
- ztot += p1.z;
- hits++;
- }
- if(hits) {
- intPt[1].z = ztot / hits;
- }
+ cerr << " q1inP && p1inQ" << endl;
#endif
+ // if pts are equal Z is chosen arbitrarily
+ intPt[0] = zGetOrInterpolateCopy(q1, p1, p2);
+ intPt[1] = zGetOrInterpolateCopy(p1, q1, q2);
+
#if GEOS_DEBUG
cerr << " intPt[0]: " << intPt[0].toString() << endl;
cerr << " intPt[1]: " << intPt[1].toString() << endl;
#endif
- return (q1 == p1) && !p1q2p2 && !q1p2q2 ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
+ return (q1 == p1) && !q2inP && !p2inQ ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
}
- if(p1q1p2 && q1p2q2) {
+ if(q1inP && p2inQ) {
#if GEOS_DEBUG
- cerr << " p1q1p2 && q1p2q2" << endl;
-#endif
- intPt[0] = q1;
-#if COMPUTE_Z
- ztot = 0;
- hits = 0;
- q1z = interpolateZ(q1, p1, p2);
- if(!std::isnan(q1z)) {
- ztot += q1z;
- hits++;
- }
- if(!std::isnan(q1.z)) {
- ztot += q1.z;
- hits++;
- }
- if(hits) {
- intPt[0].z = ztot / hits;
- }
-#endif
- intPt[1] = p2;
-#if COMPUTE_Z
- ztot = 0;
- hits = 0;
- p2z = interpolateZ(p2, q1, q2);
- if(!std::isnan(p2z)) {
- ztot += p2z;
- hits++;
- }
- if(!std::isnan(p2.z)) {
- ztot += p2.z;
- hits++;
- }
- if(hits) {
- intPt[1].z = ztot / hits;
- }
+ cerr << " q1inP && p2inQ" << endl;
#endif
+ // if pts are equal Z is chosen arbitrarily
+ intPt[0] = zGetOrInterpolateCopy(q1, p1, p2);
+ intPt[1] = zGetOrInterpolateCopy(p2, q1, q2);
+
#if GEOS_DEBUG
cerr << " intPt[0]: " << intPt[0].toString() << endl;
cerr << " intPt[1]: " << intPt[1].toString() << endl;
#endif
- return (q1 == p2) && !p1q2p2 && !q1p1q2 ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
+ return (q1 == p2) && !q2inP && !p1inQ ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
}
- if(p1q2p2 && q1p1q2) {
+ if(q2inP && p1inQ) {
#if GEOS_DEBUG
- cerr << " p1q2p2 && q1p1q2" << endl;
-#endif
- intPt[0] = q2;
-#if COMPUTE_Z
- ztot = 0;
- hits = 0;
- q2z = interpolateZ(q2, p1, p2);
- if(!std::isnan(q2z)) {
- ztot += q2z;
- hits++;
- }
- if(!std::isnan(q2.z)) {
- ztot += q2.z;
- hits++;
- }
- if(hits) {
- intPt[0].z = ztot / hits;
- }
-#endif
- intPt[1] = p1;
-#if COMPUTE_Z
- ztot = 0;
- hits = 0;
- p1z = interpolateZ(p1, q1, q2);
- if(!std::isnan(p1z)) {
- ztot += p1z;
- hits++;
- }
- if(!std::isnan(p1.z)) {
- ztot += p1.z;
- hits++;
- }
- if(hits) {
- intPt[1].z = ztot / hits;
- }
+ cerr << " q2inP && p1inQ" << endl;
#endif
+ // if pts are equal Z is chosen arbitrarily
+ intPt[0] = zGetOrInterpolateCopy(q2, p1, p2);
+ intPt[1] = zGetOrInterpolateCopy(p1, q1, q2);
#if GEOS_DEBUG
cerr << " intPt[0]: " << intPt[0].toString() << endl;
cerr << " intPt[1]: " << intPt[1].toString() << endl;
#endif
- return (q2 == p1) && !p1q1p2 && !q1p2q2 ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
+ return (q2 == p1) && !q1inP && !p2inQ ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
}
- if(p1q2p2 && q1p2q2) {
+ if(q2inP && p2inQ) {
#if GEOS_DEBUG
- cerr << " p1q2p2 && q1p2q2" << endl;
-#endif
- intPt[0] = q2;
-#if COMPUTE_Z
- ztot = 0;
- hits = 0;
- q2z = interpolateZ(q2, p1, p2);
- if(!std::isnan(q2z)) {
- ztot += q2z;
- hits++;
- }
- if(!std::isnan(q2.z)) {
- ztot += q2.z;
- hits++;
- }
- if(hits) {
- intPt[0].z = ztot / hits;
- }
-#endif
- intPt[1] = p2;
-#if COMPUTE_Z
- ztot = 0;
- hits = 0;
- p2z = interpolateZ(p2, q1, q2);
- if(!std::isnan(p2z)) {
- ztot += p2z;
- hits++;
- }
- if(!std::isnan(p2.z)) {
- ztot += p2.z;
- hits++;
- }
- if(hits) {
- intPt[1].z = ztot / hits;
- }
+ cerr << " q2inP && p2inQ" << endl;
#endif
+ // if pts are equal Z is chosen arbitrarily
+ intPt[0] = zGetOrInterpolateCopy(q2, p1, p2);
+ intPt[1] = zGetOrInterpolateCopy(p2, q1, q2);
#if GEOS_DEBUG
cerr << " intPt[0]: " << intPt[0].toString() << endl;
cerr << " intPt[1]: " << intPt[1].toString() << endl;
#endif
- return (q2 == p2) && !p1q1p2 && !q1p1q2 ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
+ return (q2 == p2) && !q1inP && !p1inQ ? POINT_INTERSECTION : COLLINEAR_INTERSECTION;
}
return NO_INTERSECTION;
}
@@ -874,24 +600,6 @@ LineIntersector::intersection(const Coordinate& p1, const Coordinate& p2,
precisionModel->makePrecise(intPtOut);
}
-
-#if COMPUTE_Z
- double ztot = 0;
- double zvals = 0;
- double zp = interpolateZ(intPtOut, p1, p2);
- double zq = interpolateZ(intPtOut, q1, q2);
- if(!std::isnan(zp)) {
- ztot += zp;
- zvals++;
- }
- if(!std::isnan(zq)) {
- ztot += zq;
- zvals++;
- }
- if(zvals) {
- intPtOut.z = ztot / zvals;
- }
-#endif // COMPUTE_Z
return intPtOut;
}
@@ -916,6 +624,142 @@ LineIntersector::intersectionSafe(const Coordinate& p1, const Coordinate& p2,
return ptInt;
}
+/* private static */
+double
+LineIntersector::zGet(const Coordinate& p, const Coordinate& q)
+{
+ double z = p.z;
+ if ( std::isnan(z) ) {
+ z = q.z; // may be NaN
+ }
+ return z;
+}
+
+/* private static */
+double
+LineIntersector::zInterpolate(const Coordinate& p, const Coordinate& p1, const Coordinate& p2)
+{
+#if GEOS_DEBUG
+ cerr << " zInterpolate " << p << ", " << p1 << ", " << p2 << std::endl;
+#endif
+ double p1z = p1.z;
+ double p2z = p2.z;
+ if (std::isnan(p1z)) {
+ return p2z; // may be NaN
+ }
+ if (std::isnan(p2z)) {
+ return p1z; // may be NaN
+ }
+ if (p.equals2D(p1)) {
+ return p1z; // not NaN
+ }
+ if (p.equals2D(p2)) {
+ return p2z; // not NaN
+ }
+ double dz = p2z - p1z;
+ if (dz == 0.0) {
+ return p1z;
+ }
+#if GEOS_DEBUG
+ cerr << " Interpolating Z from distance of p along p1-p2" << endl;
+#endif
+ // interpolate Z from distance of p along p1-p2
+ double dx = (p2.x - p1.x);
+ double dy = (p2.y - p1.y);
+ // seg has non-zero length since p1 < p < p2
+ double seglen = (dx * dx + dy * dy);
+ double xoff = (p.x - p1.x);
+ double yoff = (p.y - p1.y);
+ double plen = (xoff * xoff + yoff * yoff);
+ double frac = sqrt(plen / seglen);
+ double zoff = dz * frac;
+ double zInterpolated = p1z + zoff;
+#if GEOS_DEBUG
+ cerr << " interpolated Z: " << zInterpolated << std::endl;
+#endif
+ return zInterpolated;
+}
+
+
+/* private static */
+double
+LineIntersector::zGetOrInterpolate(const Coordinate& p, const Coordinate& p1, const Coordinate& p2)
+{
+ double z = p.z;
+ if (! std::isnan(z) ) return z;
+ return zInterpolate(p, p1, p2); // may be NaN
+}
+
+/* private static */
+Coordinate
+LineIntersector::zGetOrInterpolateCopy(
+ const Coordinate& p,
+ const Coordinate& p1,
+ const Coordinate& p2)
+{
+ Coordinate pCopy = p;
+ double z = zGetOrInterpolate(p, p1, p2);
+ pCopy.z = z;
+ return pCopy;
+}
+
+
+double
+LineIntersector::zInterpolate(const Coordinate& p,
+ const Coordinate& p1,
+ const Coordinate& p2,
+ const Coordinate& q1,
+ const Coordinate& q2)
+{
+#if GEOS_DEBUG
+ cerr << " zInterpolate(5 coords) called" << endl;
+#endif // GEOS_DEBUG
+ double zp = zInterpolate(p, p1, p2);
+#if GEOS_DEBUG
+ cerr << " zp: " << zp << endl;
+#endif // GEOS_DEBUG
+ double zq = zInterpolate(p, q1, q2);
+#if GEOS_DEBUG
+ cerr << " zq: " << zq << endl;
+#endif // GEOS_DEBUG
+ if (std::isnan(zp)) {
+ return zq; // may be NaN
+ }
+ if (std::isnan(zq)) {
+ return zp; // may be NaN
+ }
+#if GEOS_DEBUG
+ cerr << " averaging Z" << endl;
+#endif // GEOS_DEBUG
+ // both Zs have values, so average them
+ return (zp + zq) / 2.0;
+}
+
+
+/* private static */
+Coordinate
+LineIntersector::nearestEndpoint(const Coordinate& p1, const Coordinate& p2,
+ const Coordinate& q1, const Coordinate& q2)
+{
+ const Coordinate* nearestPt = &p1;
+ double minDist = Distance::pointToSegment(p1, q1, q2);
+
+ double dist = Distance::pointToSegment(p2, q1, q2);
+ if(dist < minDist) {
+ minDist = dist;
+ nearestPt = &p2;
+ }
+ dist = Distance::pointToSegment(q1, p1, p2);
+ if(dist < minDist) {
+ minDist = dist;
+ nearestPt = &q1;
+ }
+ dist = Distance::pointToSegment(q2, p1, p2);
+ if(dist < minDist) {
+ nearestPt = &q2;
+ }
+ return *nearestPt;
+}
} // namespace geos.algorithm
} // namespace geos
diff --git a/tests/unit/Makefile.am b/tests/unit/Makefile.am
index 5949522..edda8c1 100644
--- a/tests/unit/Makefile.am
+++ b/tests/unit/Makefile.am
@@ -62,6 +62,7 @@ geos_unit_SOURCES = \
algorithm/PointLocatorTest.cpp \
algorithm/RobustLineIntersectionTest.cpp \
algorithm/RobustLineIntersectorTest.cpp \
+ algorithm/RobustLineIntersectorZTest.cpp \
capi/GEOSBufferTest.cpp \
capi/GEOSBuildAreaTest.cpp \
capi/GEOSCAPIDefinesTest.cpp \
@@ -203,6 +204,7 @@ geos_unit_SOURCES = \
operation/overlayng/OverlayNGSnappingOneTest.cpp \
operation/overlayng/OverlayNGStrictModeTest.cpp \
operation/overlayng/OverlayNGTest.cpp \
+ operation/overlayng/OverlayNGZTest.cpp \
operation/overlayng/PrecisionReducerTest.cpp \
operation/overlayng/PrecisionUtilTest.cpp \
operation/overlayng/UnaryUnionNGTest.cpp \
diff --git a/tests/unit/algorithm/RobustLineIntersectorZTest.cpp b/tests/unit/algorithm/RobustLineIntersectorZTest.cpp
new file mode 100644
index 0000000..08733de
--- /dev/null
+++ b/tests/unit/algorithm/RobustLineIntersectorZTest.cpp
@@ -0,0 +1,317 @@
+//
+// Ported from JTS junit/algorithm/RobustLineIntersectorZTest.java e85881c1c179a49976fe38f0dbbff7d9fa61e65a
+
+#include <tut/tut.hpp>
+#include <utility.h>
+// geos
+#include <geos/algorithm/LineIntersector.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/LineSegment.h>
+// std
+#include <sstream>
+#include <string>
+#include <memory>
+
+
+using namespace geos::geom;
+
+typedef geos::algorithm::LineIntersector RobustLineIntersector;
+
+#define DoubleNaN std::numeric_limits<double>::quiet_NaN()
+
+namespace tut {
+//
+// Test Group
+//
+
+struct test_robustlineintersectorz_data {
+
+ void checkIntersection(const LineSegment& line1, const LineSegment& line2,
+ const Coordinate& p1, const Coordinate& p2)
+ {
+ checkIntersectionDir(line1, line2, p1, p2);
+ checkIntersectionDir(line2, line1, p1, p2);
+ LineSegment line1Rev(line1.p1, line1.p0);
+ LineSegment line2Rev(line2.p1, line2.p0);
+ checkIntersectionDir(line1Rev, line2Rev, p1, p2);
+ checkIntersectionDir(line2Rev, line1Rev, p1, p2);
+ }
+
+
+ void checkIntersectionDir(
+ const LineSegment& line1,
+ const LineSegment& line2,
+ const Coordinate& p1,
+ const Coordinate& p2)
+ {
+ RobustLineIntersector li;
+ li.computeIntersection(
+ line1.p0, line1.p1,
+ line2.p0, line2.p1);
+
+ ensure_equals(li.getIntersectionNum(), 2u);
+
+ Coordinate actual1 = li.getIntersection(0);
+ Coordinate actual2 = li.getIntersection(1);
+ // normalize actual results
+ if (actual1.compareTo(actual2) > 0) {
+ actual1 = li.getIntersection(1);
+ actual2 = li.getIntersection(0);
+ }
+
+ ensure_equals_xyz( actual1, p1 );
+ ensure_equals_xyz( actual2, p2 );
+ }
+
+ void checkIntersection(
+ const LineSegment& line1,
+ const LineSegment& line2,
+ const Coordinate& pt)
+ {
+ checkIntersectionDir(line1, line2, pt);
+ checkIntersectionDir(line2, line1, pt);
+ LineSegment line1Rev(line1.p1, line1.p0);
+ LineSegment line2Rev(line2.p1, line2.p0);
+ checkIntersectionDir(line1Rev, line2Rev, pt);
+ checkIntersectionDir(line2Rev, line1Rev, pt);
+ }
+
+ void checkIntersectionDir(
+ const LineSegment& line1,
+ const LineSegment& line2,
+ const Coordinate& pt)
+ {
+ RobustLineIntersector li;
+ li.computeIntersection(
+ line1.p0, line1.p1,
+ line2.p0, line2.p1);
+ ensure_equals(li.getIntersectionNum(), 1u);
+ Coordinate actual = li.getIntersection(0);
+ ensure_equals_xyz( actual, pt );
+ }
+
+ static Coordinate pt(double x, double y, double z) {
+ return Coordinate(x, y, z);
+ }
+
+ static Coordinate pt(double x, double y) {
+ return Coordinate(x, y);
+ }
+
+ static LineSegment line(double x1, double y1, double z1,
+ double x2, double y2, double z2)
+ {
+ return LineSegment(Coordinate(x1, y1, z1),
+ Coordinate(x2, y2, z2));
+ }
+
+ static LineSegment line(double x1, double y1,
+ double x2, double y2)
+ {
+ return LineSegment(Coordinate(x1, y1), Coordinate(x2, y2));
+ }
+
+};
+
+typedef test_group<test_robustlineintersectorz_data> group;
+typedef group::object object;
+
+group test_robustlineintersectorz_group(
+ "geos::algorithm::RobustLineIntersectorZ");
+
+
+
+
+//
+// Test Cases
+//
+
+
+// 1 - testInterior
+template<>
+template<>
+void object::test<1>
+()
+{
+ checkIntersection(
+ line(1, 1, 1, 3, 3, 3),
+ line(1, 3, 10, 3, 1, 30),
+ pt(2, 2, 11)
+ );
+}
+
+
+// 2 - testInterior2D
+template<>
+template<>
+void object::test<2>
+()
+{
+ checkIntersection(
+ line(1, 1, 3, 3),
+ line(1, 3, 3, 1),
+ pt(2, 2)
+ );
+}
+
+// testInterior3D2D
+template<>
+template<>
+void object::test<3>
+()
+{
+ checkIntersection( line(1, 1, 1, 3, 3, 3), line(1, 3, 3, 1),
+ pt(2, 2, 2));
+}
+
+// testInterior2D3D
+template<>
+template<>
+void object::test<4>
+()
+{
+ checkIntersection( line(1, 1, 3, 3), line(1, 3, 10, 3, 1, 30),
+ pt(2, 2, 20));
+}
+
+// testInterior2D3DPart
+ // result is average of line1 interpolated and line2 p0 Z
+template<>
+template<>
+void object::test<5>
+()
+{
+ checkIntersection(
+ line(1, 1, 1, 3, 3, 3),
+ line(1, 3, 10, 3, 1, DoubleNaN),
+ pt(2, 2, 6));
+}
+
+// testEndpoint
+template<>
+template<>
+void object::test<6>
+()
+{
+ checkIntersection( line(1, 1, 1, 3, 3, 3), line(3, 3, 3, 3, 1, 30),
+ pt(3, 3, 3));
+}
+
+// testEndpoint2D
+template<>
+template<>
+void object::test<7>
+()
+{
+ checkIntersection( line(1, 1, 3, 3), line(3, 3, 3, 1),
+ pt(3, 3, DoubleNaN));
+}
+
+// testEndpoint2D3D
+template<>
+template<>
+void object::test<8>
+()
+{
+ // result Z is from 3D point
+ checkIntersection( line(1, 1, 1, 3, 3, 3), line(3, 3, 3, 1),
+ pt(3, 3, 3));
+}
+
+// testInteriorEndpoint
+template<>
+template<>
+void object::test<9>
+()
+{
+ // result Z is from 3D point
+ checkIntersection( line(1, 1, 1, 3, 3, 3), line(2, 2, 10, 3, 1, 30),
+ pt(2, 2, 10));
+}
+
+// testInteriorEndpoint3D2D
+template<>
+template<>
+void object::test<10>
+()
+{
+ // result Z is interpolated
+ checkIntersection( line(1, 1, 1, 3, 3, 3), line(2, 2, 3, 1),
+ pt(2, 2, 2));
+}
+
+// testInteriorEndpoint2D3D
+ // result Z is from 3D point
+template<>
+template<>
+void object::test<11>
+()
+{
+ checkIntersection( line(1, 1, 3, 3), line(2, 2, 10, 3, 1, 20),
+ pt(2, 2, 10));
+}
+
+// testCollinearEqual
+template<>
+template<>
+void object::test<12>
+()
+{
+ checkIntersection( line(1, 1, 1, 3, 3, 3), line(1, 1, 1, 3, 3, 3),
+ pt(1, 1, 1), pt( 3, 3, 3));
+}
+
+// testCollinearEqual3D2D
+template<>
+template<>
+void object::test<13>
+()
+{
+ checkIntersection( line(1, 1, 1, 3, 3, 3), line(1, 1, 3, 3),
+ pt(1, 1, 1), pt( 3, 3, 3));
+}
+
+// testCollinearEndpoint
+template<>
+template<>
+void object::test<14>
+()
+{
+ checkIntersection( line(1, 1, 1, 3, 3, 3), line(3, 3, 3, 5, 5, 5),
+ pt(3, 3, 3));
+}
+
+// testCollinearEndpoint3D2D
+ // result Z is from 3D point
+template<>
+template<>
+void object::test<15>
+()
+{
+ checkIntersection( line(1, 1, 1, 3, 3, 3), line(3, 3, 5, 5),
+ pt(3, 3, 3));
+}
+
+// testCollinearContained
+template<>
+template<>
+void object::test<16>
+()
+{
+ checkIntersection( line(1, 1, 1, 5, 5, 5), line(3, 3, 3, 4, 4, 4),
+ pt(3, 3, 3), pt(4, 4, 4));
+}
+
+// testCollinearContained3D2D
+template<>
+template<>
+void object::test<17>
+()
+{
+ // result Z is interpolated
+ checkIntersection( line(1, 1, 1, 5, 5, 5), line(3, 3, 4, 4),
+ pt(3, 3, 3), pt(4, 4, 4));
+}
+
+} // namespace tut
+
diff --git a/tests/unit/operation/overlayng/OverlayNGZTest.cpp b/tests/unit/operation/overlayng/OverlayNGZTest.cpp
new file mode 100644
index 0000000..6ee56af
--- /dev/null
+++ b/tests/unit/operation/overlayng/OverlayNGZTest.cpp
@@ -0,0 +1,56 @@
+//
+// Test Suite for handling of Z in geos::operation::overlayng::OverlayNG class.
+
+#include <tut/tut.hpp>
+#include <utility.h>
+
+// geos
+#include <geos/operation/overlayng/OverlayNG.h>
+
+// std
+#include <memory>
+
+using namespace geos::geom;
+using namespace geos::operation::overlayng;
+using geos::io::WKTReader;
+
+namespace tut {
+//
+// Test Group
+//
+
+// Common data used by all tests
+struct test_overlayngz_data {
+
+ WKTReader r;
+
+ void checkIntersection(const std::string& wktA, const std::string& wktB, const std::string& wktExpected)
+ {
+ std::unique_ptr<Geometry> a = r.read(wktA);
+ std::unique_ptr<Geometry> b = r.read(wktB);
+ std::unique_ptr<Geometry> expected = r.read(wktExpected);
+ std::unique_ptr<Geometry> result = OverlayNG::overlay(a.get(), b.get(), OverlayNG::INTERSECTION);
+ ensure_equals_geometry(expected.get(), result.get());
+ }
+
+};
+
+typedef test_group<test_overlayngz_data> group;
+typedef group::object object;
+
+group test_overlayngz_group("geos::operation::overlayng::OverlayNGZ");
+
+//
+// Test Cases
+//
+
+// testLineIntersectionInterpolated
+template<>
+template<>
+void object::test<1> ()
+{
+ checkIntersection("LINESTRING (0 0 0, 10 10 10)", "LINESTRING (10 0 0, 0 10 10)", "POINT(5 5 5)");
+}
+
+
+} // namespace tut
diff --git a/tests/unit/utility.h b/tests/unit/utility.h
index 40db8bc..c3da601 100644
--- a/tests/unit/utility.h
+++ b/tests/unit/utility.h
@@ -92,6 +92,21 @@ instanceOf(InstanceType const* instance)
return dynamic_cast<Type const*>(instance);
}
+inline void
+ensure_equals_xyz(geos::geom::Coordinate const& actual,
+ geos::geom::Coordinate const& expected)
+{
+ ensure_equals("Coordinate X", actual.x, expected.x );
+ ensure_equals("Coordinate Y", actual.y, expected.y );
+ if ( std::isnan(expected.z) ) {
+ ensure("Coordinate Z should be NaN", std::isnan(actual.z) );
+ } else {
+ ensure_equals("Coordinate Z", actual.z, expected.z );
+ }
+}
+
+
+
//
// Geometries structure comparators
//
-----------------------------------------------------------------------
Summary of changes:
include/geos/algorithm/LineIntersector.h | 48 ++
src/algorithm/LineIntersector.cpp | 596 ++++++++-------------
tests/unit/Makefile.am | 2 +
.../unit/algorithm/RobustLineIntersectorZTest.cpp | 317 +++++++++++
tests/unit/operation/overlayng/OverlayNGZTest.cpp | 56 ++
tests/unit/utility.h | 15 +
6 files changed, 658 insertions(+), 376 deletions(-)
create mode 100644 tests/unit/algorithm/RobustLineIntersectorZTest.cpp
create mode 100644 tests/unit/operation/overlayng/OverlayNGZTest.cpp
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list