[geos-commits] [SCM] GEOS branch main updated. 11f1851a82bdf5f9983f03d1b967adfb5da6fb39
git at osgeo.org
git at osgeo.org
Tue Jun 16 11:27:15 PDT 2026
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 11f1851a82bdf5f9983f03d1b967adfb5da6fb39 (commit)
via 47a2a8fba3190aa8604f2e2ef55e649cd67f3f7a (commit)
via 420c2a988d54fc07e9ff8a513f7df1fdca7ad741 (commit)
via 3da292784f83e0bfce769442d4d7e14cdc964bc6 (commit)
via f0ffac4dc49d012f3d80cb28f76e9ff10d79772f (commit)
via fac5376f75356ad9cd8b6f08ef0336f1d56e006c (commit)
via c06102e28158cd1393a68c767f36087b249e4ee0 (commit)
via f14c453ade3f86a98b3dafeaf5cca750593d4e70 (commit)
via dc3935ba35f5fd6328bddb3637c7e0eeba307dd8 (commit)
via 1aeef7493bc6d0804d094021fc1013a7227791af (commit)
via b0506ea14bdc0339fb8ab61eb95ba5ec26762332 (commit)
via 43edb6d35a6c66fe82712033bb3abd53f89fce3a (commit)
via 095c72709ea5043b4a1433d2adb3611410564cd3 (commit)
via fa14b9a717e2231679859d20adae502ebff07492 (commit)
from ac5da1cb2ed3e76cd286bb0f715e7fccde89281e (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 11f1851a82bdf5f9983f03d1b967adfb5da6fb39
Author: Daniel Baston <dbaston at gmail.com>
Date: Mon Jun 8 08:49:21 2026 -0400
Apply suggestions from @rouault
diff --git a/include/geos/operation/overlayng/Edge.h b/include/geos/operation/overlayng/Edge.h
index b88a44fe7..3ee964c5d 100644
--- a/include/geos/operation/overlayng/Edge.h
+++ b/include/geos/operation/overlayng/Edge.h
@@ -64,7 +64,7 @@ private:
int bDepthDelta = 0;
bool bIsHole = false;
std::shared_ptr<const geom::CoordinateSequence> pts;
- bool ptsCurved;
+ const bool ptsCurved;
// Methods
diff --git a/src/operation/overlayng/OverlayEdgeRing.cpp b/src/operation/overlayng/OverlayEdgeRing.cpp
index 12f4c8e92..98db9f536 100644
--- a/src/operation/overlayng/OverlayEdgeRing.cpp
+++ b/src/operation/overlayng/OverlayEdgeRing.cpp
@@ -284,7 +284,7 @@ OverlayEdgeRing::isPointInOrOut(const OverlayEdgeRing& otherRing) const {
struct PointTester : public geom::CoordinateFilter {
public:
- explicit PointTester(const OverlayEdgeRing* p_ring) : m_er(p_ring), m_result(false) {}
+ explicit PointTester(const OverlayEdgeRing* p_ring) : m_er(p_ring) {}
void filter_ro(const CoordinateXY* pt) override {
Location loc = m_er->locate(*pt);
@@ -311,7 +311,7 @@ OverlayEdgeRing::isPointInOrOut(const OverlayEdgeRing& otherRing) const {
private:
const OverlayEdgeRing* m_er;
bool m_done{false};
- bool m_result;
+ bool m_result{false};
};
PointTester tester(this);
commit 47a2a8fba3190aa8604f2e2ef55e649cd67f3f7a
Author: Daniel Baston <dbaston at gmail.com>
Date: Sun Jun 7 14:42:06 2026 -0400
CircularArcIntersector: exactly preserve endpoints in cocircular arc intersection
diff --git a/src/algorithm/CircularArcIntersector.cpp b/src/algorithm/CircularArcIntersector.cpp
index 8d3fabf58..5f82d3230 100644
--- a/src/algorithm/CircularArcIntersector.cpp
+++ b/src/algorithm/CircularArcIntersector.cpp
@@ -427,32 +427,34 @@ CircularArcIntersector::addCocircularIntersection(double startAngle, double endA
CoordinateXYZM computedMidPt(CircularArcs::createPoint(center, radius, theta1));
CoordinateXYZM computedEndPt(CircularArcs::createPoint(center, radius, endAngle));
- if (precisionModel) {
- precisionModel->makePrecise(computedStartPt);
- precisionModel->makePrecise(computedMidPt);
- precisionModel->makePrecise(computedEndPt);
- }
-
// Check to see if the endpoints of the intersection match the endpoints of either of
// the input arcs. Use angles for the check to avoid missing an endpoint intersection from
// inaccuracy in the point construction.
- if (startAngle == arc1.theta0()) {
+ if (startAngle == Angle::normalizePositive(arc1.theta0())) {
+ computedStartPt = arc1.p0();
setFromEndpoint(computedStartPt, arc1, 0);
- } else if (startAngle == arc1.theta2()) {
+ } else if (startAngle == Angle::normalizePositive(arc1.theta2())) {
+ computedStartPt = arc1.p2();
setFromEndpoint(computedStartPt, arc1, 2);
- } else if (startAngle == arc2.theta0()) {
+ } else if (startAngle == Angle::normalizePositive(arc2.theta0())) {
+ computedStartPt = arc2.p0();
setFromEndpoint(computedStartPt, arc2, 0);
- } else if (startAngle == arc2.theta2()) {
+ } else if (startAngle == Angle::normalizePositive(arc2.theta2())) {
+ computedStartPt = arc2.p2();
setFromEndpoint(computedStartPt, arc2, 2);
}
- if (endAngle == arc1.theta0()) {
+ if (endAngle == Angle::normalizePositive(arc1.theta0())) {
+ computedEndPt = arc1.p0();
setFromEndpoint(computedEndPt, arc1, 0);
- } else if (endAngle == arc1.theta2()) {
+ } else if (endAngle == Angle::normalizePositive(arc1.theta2())) {
+ computedEndPt = arc1.p2();
setFromEndpoint(computedEndPt, arc1, 2);
- } else if (endAngle == arc2.theta0()) {
+ } else if (endAngle == Angle::normalizePositive(arc2.theta0())) {
+ computedEndPt = arc2.p0();
setFromEndpoint(computedEndPt, arc2, 0);
- } else if (endAngle == arc2.theta2()) {
+ } else if (endAngle == Angle::normalizePositive(arc2.theta2())) {
+ computedEndPt = arc2.p2();
setFromEndpoint(computedEndPt, arc2, 2);
}
@@ -460,6 +462,12 @@ CircularArcIntersector::addCocircularIntersection(double startAngle, double endA
interpolateZM(arc1, arc2, computedMidPt);
interpolateZM(arc1, arc2, computedEndPt);
+ if (precisionModel) {
+ precisionModel->makePrecise(computedStartPt);
+ precisionModel->makePrecise(computedMidPt);
+ precisionModel->makePrecise(computedEndPt);
+ }
+
auto seq = std::make_unique<CoordinateSequence>(3, constructZ, constructM);
seq->setAt(computedStartPt, 0);
seq->setAt(computedMidPt, 1);
diff --git a/tests/unit/algorithm/CircularArcIntersectorTest.cpp b/tests/unit/algorithm/CircularArcIntersectorTest.cpp
index e9899213d..5d0123c17 100644
--- a/tests/unit/algorithm/CircularArcIntersectorTest.cpp
+++ b/tests/unit/algorithm/CircularArcIntersectorTest.cpp
@@ -1733,6 +1733,26 @@ void object::test<86>()
}
+template<>
+template<>
+void object::test<87>()
+{
+ set_test_name("arc intersected with itself yields arc with identical endpoints");
+
+ CircularArcIntersector cai;
+
+ auto arc = CircularArc::create(XY{0, 0}, XY{2, 0}, XY{2, 1});
+
+ cai.intersects(arc, arc);
+
+ ensure_equals(cai.getNumArcs(), 1u);
+
+ const CircularArc& arcOut = cai.getArc(0);
+
+ ensure_equals(arc.p0(), arcOut.p0());
+ ensure_equals(arc.p2(), arcOut.p2());
+}
+
// TODO: check Z values of arc result centerpoints
// TODO: add tests for seg/seg
commit 420c2a988d54fc07e9ff8a513f7df1fdca7ad741
Author: Daniel Baston <dbaston at gmail.com>
Date: Mon Jun 8 08:35:41 2026 -0400
NodableArcString: properly handle splits at arc endpoints
diff --git a/src/noding/NodableArcString.cpp b/src/noding/NodableArcString.cpp
index 4b824b3fa..d05865e79 100644
--- a/src/noding/NodableArcString.cpp
+++ b/src/noding/NodableArcString.cpp
@@ -21,6 +21,8 @@
using geos::geom::CoordinateXYZM;
using geos::geom::CircularArc;
+#define DEBUG_NODABLE_ARC_STRING 0
+
namespace geos::noding {
static double
@@ -41,12 +43,22 @@ NodableArcString::NodableArcString(std::vector<geom::CircularArc> arcs, const st
{
}
-static
-std::vector<CoordinateXYZM> prepareArcPoints(const CircularArc& arc, std::vector<CoordinateXYZM> splitPoints)
+struct SplitPoints {
+
+ std::vector<CoordinateXYZM> points{};
+ bool splitStart{false};
+ bool splitEnd{false};
+};
+
+static SplitPoints
+prepareArcPoints(const CircularArc& arc, std::vector<CoordinateXYZM> splitPoints)
{
const bool isCCW = arc.getOrientation() == algorithm::Orientation::COUNTERCLOCKWISE;
const geom::CoordinateXY& center = arc.getCenter();
+ bool splitStart = false;
+ bool splitEnd = false;
+
// Some potential split points may be skipped, for example, if they would create an arc section that is
// too short to have a constructed midpoint. Because the results of this logic could depend on the
// direction in which the arc is processed, we reverse clockwise arcs and then reverse the list of
@@ -65,7 +77,17 @@ std::vector<CoordinateXYZM> prepareArcPoints(const CircularArc& arc, std::vector
double pa0 = geom::Quadrant::pseudoAngle(center, p0);
double pa1 = geom::Quadrant::pseudoAngle(center, p1);
- return pseudoAngleDiffCCW(paStart, pa0) < pseudoAngleDiffCCW(paStart, pa1);
+ double diff0 = pseudoAngleDiffCCW(paStart, pa0);
+ double diff1 = pseudoAngleDiffCCW(paStart, pa1);
+
+ if (diff0 < diff1) {
+ return true;
+ }
+ if (diff0 > diff1) {
+ return false;
+ }
+
+ return p0 < p1;
});
// Add ending point of input arc
@@ -75,20 +97,50 @@ std::vector<CoordinateXYZM> prepareArcPoints(const CircularArc& arc, std::vector
splitPoints.push_back(p2);
}
+#if DEBUG_NODABLE_ARC_STRING
+ std::cout << std::setprecision(17);
+ std::cout << std::endl;
+ std::cout << "Edge from " << arc.p0() << " to " << arc.p2() << " has " << splitPoints.size() << " potential split points:" << std::endl;
+ std::cout << "Edge is CIRCULARSTRING " << *arc.getCoordinateSequence() << std::endl;
+ std::cout << "Start point is " << retained.back() << " pa " << paStart << std::endl;
+ std::cout << "Potential split points:" << std::endl;
+ for (std::size_t i = 0; i < splitPoints.size(); i++) {
+ std::cout << " " << splitPoints[i] << " paDiff " << pseudoAngleDiffCCW(paStart, geom::Quadrant::pseudoAngle(center, splitPoints[i]));
+ if (i > 0) {
+ std::cout << " " << splitPoints[i].distance(splitPoints[i-1]);
+ }
+ std::cout << std::endl;
+ }
+#endif
+
for (const auto& p2 : splitPoints) {
auto& p0 = retained.back();
+#if DEBUG_NODABLE_ARC_STRING
+ std::cout << "p0 = " << p0 << std::endl;
+#endif
+
if (p2.equals2D(p0)) {
+#if DEBUG_NODABLE_ARC_STRING
+ std::cout << "Split point " << p2 << " equal to start point " << p0 << std::endl;
+#endif
if (std::isnan(p0.z) && !std::isnan(p2.z)) {
p0.z = p2.z;
}
if (std::isnan(p0.m) && !std::isnan(p2.m)) {
p0.m = p2.m;
}
+
+ if (retained.size() == 1) {
+ splitStart = true;
+ }
continue;
}
if (!arc.containsPointOnCircle(p2)) {
+#if DEBUG_NODABLE_ARC_STRING
+ std::cout << "Skipping split point " << p2 << " because it is outside of the arc" << std::endl;
+#endif
continue;
}
@@ -98,6 +150,22 @@ std::vector<CoordinateXYZM> prepareArcPoints(const CircularArc& arc, std::vector
const geom::CoordinateXY p1 = algorithm::CircularArcs::getMidpoint(p0, p2, center, arc.getRadius(), true);
if (p1.equals2D(p0) || p1.equals2D(p2)) {
+#if DEBUG_NODABLE_ARC_STRING
+ std::cout << "Skipping split point " << p2 << " because the calculated arc midpoint " << p1 << " equals one of the endpoints" << std::endl;
+#endif
+ if (retained.size() == 1) {
+ splitStart = true;
+ }
+ continue;
+ }
+
+ if (algorithm::Orientation::index(p1, p0, p2) == algorithm::Orientation::COLLINEAR) {
+#if DEBUG_NODABLE_ARC_STRING
+ std::cout << "Skipping split point " << p2 << " because the calculated arc midpoint " << p1 << " is collinear with the endpoints" << std::endl;
+#endif
+ if (retained.size() == 1) {
+ splitStart = true;
+ }
continue;
}
@@ -106,10 +174,19 @@ std::vector<CoordinateXYZM> prepareArcPoints(const CircularArc& arc, std::vector
const double t1 = algorithm::Angle::normalizePositive(algorithm::CircularArcs::getAngle(p1, center));
const double t2 = algorithm::Angle::normalizePositive(isCCW ? arc.theta2() : arc.theta0());
- if (!algorithm::Angle::isWithinCCW(t1, t0, t2)) {
+ if (!algorithm::Angle::isWithinCCW(t1, t0, t2)) { // != isCCW) {
+#if DEBUG_NODABLE_ARC_STRING
+ std::cout << "Skipping split point " << p2 << " because the calculated arc midpoint " << p1 << " does not fall within the arc from " << p0 << " to " << p2 << std::endl;
+#endif
+ if (retained.size() == 1) {
+ splitStart = true;
+ }
continue;
}
+#if DEBUG_NODABLE_ARC_STRING
+ std::cout << "Keeping split point " << p2 << std::endl;
+#endif
retained.push_back(p2);
}
@@ -120,6 +197,7 @@ std::vector<CoordinateXYZM> prepareArcPoints(const CircularArc& arc, std::vector
CoordinateXYZM& back = retained.back();
if (!back.equals2D(p2)) {
+ splitEnd = true;
back.x = p2.x;
back.y = p2.y;
}
@@ -133,15 +211,75 @@ std::vector<CoordinateXYZM> prepareArcPoints(const CircularArc& arc, std::vector
if (!isCCW) {
std::reverse(retained.begin(), retained.end());
+ std::swap(splitStart, splitEnd);
}
- return retained;
+#if DEBUG_NODABLE_ARC_STRING
+ std::cout << "Retained split points:" << std::endl;
+ for (std::size_t i = 0; i < retained.size(); i++) {
+ std::cout << " " << retained[i] << " paDiff " << pseudoAngleDiffCCW(paStart, geom::Quadrant::pseudoAngle(center, retained[i]));
+ if (i > 0) {
+ std::cout << " " << retained[i].distance(retained[i-1]);
+ }
+ std::cout << std::endl;
+ }
+#endif
+
+ return SplitPoints{std::move(retained), splitStart, splitEnd};
}
+class ArcBuilder {
+
+public:
+
+ ArcBuilder(std::vector<std::unique_ptr<ArcString>>& splitArcs, const void* data, bool constructZ, bool constructM) :
+ m_data(data),
+ m_splitArcs(splitArcs),
+ m_constructZ(constructZ),
+ m_constructM(constructM) {}
+
+ void addArc(const CoordinateXYZM& p0, const CoordinateXYZM& p1, const CoordinateXYZM& p2, const geom::CoordinateXY& center, double radius, int orientation) {
+ if (!m_coords) {
+ m_coords = std::make_unique<geom::CoordinateSequence>(0, m_constructZ, m_constructM);
+ }
+
+ m_coords->add(p0, false);
+ m_coords->add(p1, false);
+ m_coords->add(p2, false);
+
+ m_arcs.emplace_back(*m_coords, m_coords->getSize() - 3, center, radius, orientation);
+ }
+
+ void addArc(const geom::CoordinateSequence& seq, std::size_t pos, const geom::CoordinateXY& center, double radius, int orientation) {
+ if (!m_coords) {
+ m_coords = std::make_unique<geom::CoordinateSequence>(0, m_constructZ, m_constructM);
+ }
+
+ m_coords->add(seq, pos, pos + 2, false);
+ m_arcs.emplace_back(*m_coords, m_coords->getSize() - 3, center, radius, orientation);
+ }
+
+ void finish() {
+ if (m_arcs.empty()) {
+ return;
+ }
+
+ m_splitArcs.push_back(std::make_unique<NodableArcString>(std::move(m_arcs), std::move(m_coords), m_constructZ, m_constructM, m_data));
+ m_arcs = std::vector<CircularArc>();
+ }
+
+private:
+ const void* m_data;
+ std::vector<std::unique_ptr<ArcString>>& m_splitArcs;
+ std::vector<CircularArc> m_arcs;
+ std::unique_ptr<geom::CoordinateSequence> m_coords;
+ const bool m_constructZ;
+ const bool m_constructM;
+};
+
void
NodableArcString::getNoded(std::vector<std::unique_ptr<ArcString>>& splitArcs) {
- auto dstSeq = std::make_unique<geom::CoordinateSequence>(0, m_constructZ, m_constructM);
- std::vector<geom::CircularArc> arcs;
+ ArcBuilder builder(splitArcs, getData(), m_constructZ, m_constructM);
for (size_t arcIndex = 0; arcIndex < m_arcs.size(); arcIndex++) {
const CircularArc& toSplit = m_arcs[arcIndex];
@@ -149,36 +287,28 @@ NodableArcString::getNoded(std::vector<std::unique_ptr<ArcString>>& splitArcs) {
const double radius = toSplit.getRadius();
const int orientation = toSplit.getOrientation();
- bool arcIsSplit = true;
- bool createArcString = true;
- const bool preserveControlPoint = true;
- std::vector<CoordinateXYZM> arcPoints;
- const auto it = m_adds.find(arcIndex);
- if (it == m_adds.end()) {
- arcIsSplit = false;
- createArcString = false;
- } else {
- arcPoints = prepareArcPoints(toSplit, it->second);
+ SplitPoints split;
- if (arcPoints.size() == 2) {
- // All added nodes collapsed
- // Still need to know if arc was split.
- arcIsSplit = false;
- }
+ const bool preserveControlPoint = true;
+
+ const auto it = m_adds.find(arcIndex);
+ if (it != m_adds.end()) {
+ split = prepareArcPoints(toSplit, it->second);
}
- if (preserveControlPoint && !arcIsSplit) {
+ if (split.splitStart) {
+ builder.finish();
+ }
+
+ if (preserveControlPoint && split.points.size() <= 2) {
// No nodes added, just copy the coordinates into the sequence.
const geom::CoordinateSequence* srcSeq = m_arcs[arcIndex].getCoordinateSequence();
std::size_t srcPos = m_arcs[arcIndex].getCoordinatePosition();
- dstSeq->add(*srcSeq, srcPos, srcPos + 2, false);
- std::size_t dstPos = dstSeq->getSize() - 3;
- arcs.emplace_back(*dstSeq, dstPos, center, radius, orientation);
- if (createArcString) {
- splitArcs.push_back(std::make_unique<NodableArcString>(std::move(arcs), std::move(dstSeq), m_constructZ, m_constructM, getData()));
- dstSeq = std::make_unique<geom::CoordinateSequence>(0, m_constructZ, m_constructM);
- arcs.clear();
+ builder.addArc(*srcSeq, srcPos, center, radius, orientation);
+
+ if (split.splitEnd) {
+ builder.finish();
}
continue;
@@ -186,9 +316,9 @@ NodableArcString::getNoded(std::vector<std::unique_ptr<ArcString>>& splitArcs) {
const bool isCCW = orientation == algorithm::Orientation::COUNTERCLOCKWISE;
- for (std::size_t i = 1; i < arcPoints.size(); i++) {
- const CoordinateXYZM& p0 = arcPoints[i - 1];
- const CoordinateXYZM& p2 = arcPoints[i];
+ for (std::size_t i = 1; i < split.points.size(); i++) {
+ const CoordinateXYZM& p0 = split.points[i - 1];
+ const CoordinateXYZM& p2 = split.points[i];
// TODO: Check if control point of original arc falls into this section,
// and use it instead of calculating a midpoint here?
@@ -196,28 +326,18 @@ NodableArcString::getNoded(std::vector<std::unique_ptr<ArcString>>& splitArcs) {
p1.z = (p0.z + p2.z) / 2;
p1.m = (p0.m + p2.m) / 2;
- if (dstSeq->isEmpty()) {
- dstSeq->add(p0);
- }
- dstSeq->add(p1);
- dstSeq->add(p2);
-
- const std::size_t dstPos = dstSeq->getSize() - 3;
- arcs.emplace_back(*dstSeq, dstPos, center, radius, orientation);
+ builder.addArc(p0, p1, p2, center, radius, orientation);
// Finish the ArcString, start a new one.
- const bool isSplitPoint = i != arcPoints.size() - 1;
+ const bool isSplitPoint = (i < split.points.size() - 1) || split.splitEnd;
if (isSplitPoint) {
- splitArcs.push_back(std::make_unique<NodableArcString>(std::move(arcs), std::move(dstSeq), m_constructZ, m_constructM, getData()));
- dstSeq = std::make_unique<geom::CoordinateSequence>(0, m_constructZ, m_constructM);
- arcs.clear();
+ builder.finish();
+
}
}
}
- if (!arcs.empty()) {
- splitArcs.push_back(std::make_unique<NodableArcString>(std::move(arcs), std::move(dstSeq), m_constructZ, m_constructM, getData()));
- }
+ builder.finish();
}
}
\ No newline at end of file
diff --git a/tests/unit/operation/overlayng/OverlayNGTest.cpp b/tests/unit/operation/overlayng/OverlayNGTest.cpp
index 46b8155fa..2da3cbcca 100644
--- a/tests/unit/operation/overlayng/OverlayNGTest.cpp
+++ b/tests/unit/operation/overlayng/OverlayNGTest.cpp
@@ -6,9 +6,6 @@
// geos
#include <geos/operation/overlayng/OverlayNG.h>
-#include <geos/noding/SimpleNoder.h>
-#include <geos/algorithm/CircularArcIntersector.h>
-#include <geos/noding/ArcIntersectionAdder.h>
// std
#include <memory>
@@ -818,4 +815,32 @@ void object::test<62>()
testOverlay(a, a, exp, OverlayNG::UNION, 0);
}
+template<>
+template<>
+void object::test<63>()
+{
+ set_test_name("Union of CurvePolygon / CurvePolygon -> CurvePolygon");
+
+ std::string a = "CURVEPOLYGON (COMPOUNDCURVE (CIRCULARSTRING (-1.1087328289561522 -0.2365594619955851, -1.2116499629412223 -0.2219284209572556, -1.3064159292035398 -0.1792035398230087), (-1.3064159292035398 -0.1792035398230087, -1.2982128297312812 -0.0461310372730352), CIRCULARSTRING (-1.2982128297312812 -0.0461310372730352, -1.2213760168089367 -0.1591592707276735, -1.1087328289561522 -0.2365594619955851)))";
+ std::string b = "CURVEPOLYGON (COMPOUNDCURVE ((-1.277492733105557 0.2899949746553787, -0.8433316896087282 -0.1241204385865409), CIRCULARSTRING (-0.8433316896087282 -0.1241204385865409, -0.9646131918909476 -0.2072935239179199, -1.1087328289561522 -0.2365594619955851, -1.2213760168089367 -0.1591592707276735, -1.2982128297312812 -0.0461310372730352), (-1.2982128297312812 -0.0461310372730352, -1.277492733105557 0.2899949746553787)))";
+
+ std::string exp = "CURVEPOLYGON (COMPOUNDCURVE (CIRCULARSTRING (-1.1087328289561522 -0.2365594619955851, -1.2116499629412223 -0.2219284209572556, -1.3064159292035398 -0.1792035398230087), (-1.3064159292035398 -0.1792035398230087, -1.2982128297312812 -0.0461310372730352, -1.277492733105557 0.2899949746553787, -0.8433316896087282 -0.1241204385865409), CIRCULARSTRING (-0.8433316896087282 -0.1241204385865409, -0.9646131918909476 -0.2072935239179199, -1.1087328289561522 -0.2365594619955851)))";
+
+ testOverlay(a, b, exp, OverlayNG::UNION, 0);
+}
+
+template<>
+template<>
+void object::test<64>()
+{
+ set_test_name("Union of CurvePolygon / CurvePolygon -> CurvePolygon (2)");
+
+ std::string a = "MULTISURFACE (CURVEPOLYGON (COMPOUNDCURVE (CIRCULARSTRING (-1.266592920353982 0.4668141592920354, -1.1143440097781072 0.502208020974126, -0.9610873502585282 0.4714695673823874, -1.1432212469102676 0.4224569115784076, -1.277492733105557 0.2899949746553787), (-1.277492733105557 0.2899949746553787, -1.266592920353982 0.4668141592920354))), CURVEPOLYGON (COMPOUNDCURVE (CIRCULARSTRING (-0.7885401828539679 0.3170953789649493, -0.7421261955995471 0.0873215546859087, -0.8433316896087282 -0.1241204385865409), (-0.8433316896087282 -0.1241204385865409, -0.7885401828539679 0.3170953789649493))), CURVEPOLYGON (COMPOUNDCURVE (CIRCULARSTRING (-1.1087328289561522 -0.2365594619955851, -1.2116499629412223 -0.2219284209572556, -1.3064159292035398 -0.1792035398230087), (-1.3064159292035398 -0.1792035398230087, -1.2982128297312812 -0.0461310372730352), CIRCULARSTRING (-1.2982128297312812 -0.0461310372730352, -1.2213760168089367 -0.1591592707276735, -1.1087328289561522 -0.23655946199
55851))))";
+ std::string b = "CURVEPOLYGON (COMPOUNDCURVE (CIRCULARSTRING (-0.9610873502585282 0.4714695673823874, -0.8624070192245114 0.4081497467816552, -0.7885401828539679 0.3170953789649493), (-0.7885401828539679 0.3170953789649493, -0.8433316896087282 -0.1241204385865409), CIRCULARSTRING (-0.8433316896087282 -0.1241204385865409, -0.9646131918909476 -0.2072935239179199, -1.1087328289561522 -0.2365594619955851, -1.2213760168089367 -0.1591592707276735, -1.2982128297312812 -0.0461310372730352), (-1.2982128297312812 -0.0461310372730352, -1.277492733105557 0.2899949746553787), CIRCULARSTRING (-1.277492733105557 0.2899949746553787, -1.1432212469102676 0.4224569115784076, -0.9610873502585282 0.4714695673823874)))";
+
+ std::string exp = "CURVEPOLYGON (COMPOUNDCURVE (CIRCULARSTRING (-1.266592920353982 0.4668141592920354, -1.1143440097781072 0.502208020974126, -0.9610873502585282 0.4714695673823874, -0.8624070192245114 0.4081497467816552, -0.7885401828539679 0.3170953789649493, -0.7421261955995471 0.0873215546859087, -0.8433316896087282 -0.1241204385865409, -0.9646131918909476 -0.2072935239179199, -1.1087328289561522 -0.2365594619955851, -1.2116499629412223 -0.2219284209572556, -1.3064159292035398 -0.1792035398230087), (-1.3064159292035398 -0.1792035398230087, -1.2982128297312812 -0.0461310372730352, -1.277492733105557 0.2899949746553787, -1.266592920353982 0.4668141592920354)))";
+
+ testOverlay(a, b, exp, OverlayNG::UNION, 0);
+}
+
} // namespace tut
commit 3da292784f83e0bfce769442d4d7e14cdc964bc6
Author: Daniel Baston <dbaston at gmail.com>
Date: Sun Jun 7 12:03:15 2026 -0400
CircularArcIntersector: improve robustness of arc/arc endpoint intersection
diff --git a/src/algorithm/CircularArcIntersector.cpp b/src/algorithm/CircularArcIntersector.cpp
index ad79ee190..8d3fabf58 100644
--- a/src/algorithm/CircularArcIntersector.cpp
+++ b/src/algorithm/CircularArcIntersector.cpp
@@ -235,55 +235,39 @@ CircularArcIntersector::intersects(const CircularArc& arc1, const CircularArc& a
if (a == 0 || (d == 0 && r1 == r2)) {
computeCocircularIntersection(arc1, arc2);
} else {
- // Explicitly add endpoint intersections that may be missed or inexactly computed.
- if (arc1.p0().equals2D(arc2.p0()) && !hasIntersection(arc1.p0())) {
- addArcArcIntersectionPoint(arc1.p0(), arc1, arc2);
+ // Compute interior intersection points.
+ const double dx = c2.x-c1.x;
+ const double dy = c2.y-c1.y;
+
+ // point where a line between the two circle center points intersects
+ // the radical line
+ CoordinateXY p{c1.x + a* dx/d, c1.y+a* dy/d};
+
+ // distance from p to the intersection points
+ const double h = std::sqrt(r1*r1 - a*a);
+
+ CoordinateXY isect0{p.x + h* dy/d, p.y - h* dx/d };
+ CoordinateXY isect1{p.x - h* dy/d, p.y + h* dx/d };
+
+ // Check to see if computed intersection points are inexact versions of an endpoint intersection
+ const CoordinateXY& ap0 = arc1.p0();
+ const CoordinateXY& ap2 = arc1.p2();
+ const CoordinateXY& bp0 = arc2.p0();
+ const CoordinateXY& bp2 = arc2.p2();
+
+ if (ap0 == bp0 || ap0 == bp2) {
+ closestPoint(isect0, isect1, 2, ap0) = ap0;
}
- if (arc1.p0().equals2D(arc2.p2()) && !hasIntersection(arc1.p0())) {
- addArcArcIntersectionPoint(arc1.p0(), arc1, arc2);
- }
- if (arc1.p2().equals2D(arc2.p0()) && !hasIntersection(arc1.p2())) {
- addArcArcIntersectionPoint(arc1.p2(), arc1, arc2);
- }
- if (arc1.p2().equals2D(arc2.p2()) && !hasIntersection(arc1.p2())) {
- addArcArcIntersectionPoint(arc1.p2(), arc1, arc2);
+ if (ap2 == bp0 || ap2 == bp2) {
+ closestPoint(isect0, isect1, 2, ap2) = ap2;
}
- if (nPt < 2) {
- // Compute interior intersection points.
- const double dx = c2.x-c1.x;
- const double dy = c2.y-c1.y;
+ if (arc1.containsPointOnCircle(isect0) && arc2.containsPointOnCircle(isect0)) {
+ addArcArcIntersectionPoint(isect0, arc1, arc2);
+ }
- // point where a line between the two circle center points intersects
- // the radical line
- CoordinateXY p{c1.x + a* dx/d, c1.y+a* dy/d};
-
- // distance from p to the intersection points
- const double h = std::sqrt(r1*r1 - a*a);
-
- CoordinateXY isect0{p.x + h* dy/d, p.y - h* dx/d };
- CoordinateXY isect1{p.x - h* dy/d, p.y + h* dx/d };
-
- // One of the computed intersection points may be an inexact version of an endpoint.
- // If we already have an endpoint intersection, we need to process the farther-away
- // computed point first.
- if (nPt == 1 && intPt[0].distance(isect0) < intPt[0].distance(isect1)) {
- std::swap(isect0, isect1);
- }
-
- for (const CoordinateXY& computedIntPt : {isect0, isect1}) {
- if (nPt > 0 && computedIntPt.equals2D(intPt[0])) {
- continue;
- }
-
- if (nPt > 1) {
- continue;
- }
-
- if (arc1.containsPointOnCircle(computedIntPt) && arc2.containsPointOnCircle(computedIntPt)) {
- addArcArcIntersectionPoint(computedIntPt, arc1, arc2);
- }
- }
+ if (isect1 != isect0 && arc1.containsPointOnCircle(isect1) && arc2.containsPointOnCircle(isect1)) {
+ addArcArcIntersectionPoint(isect1, arc1, arc2);
}
}
commit f0ffac4dc49d012f3d80cb28f76e9ff10d79772f
Author: Daniel Baston <dbaston at gmail.com>
Date: Sun Jun 7 11:50:32 2026 -0400
CircularArc: Add getReverseDirectionPoint
diff --git a/include/geos/geom/CircularArc.h b/include/geos/geom/CircularArc.h
index 9f9bbc693..800eb0144 100644
--- a/include/geos/geom/CircularArc.h
+++ b/include/geos/geom/CircularArc.h
@@ -110,6 +110,8 @@ public:
/// line tangent to the arc at the origin point.
CoordinateXY getDirectionPoint() const;
+ CoordinateXY getReverseDirectionPoint() const;
+
Envelope getEnvelope() const;
/// Return the length of the arc
diff --git a/src/geom/CircularArc.cpp b/src/geom/CircularArc.cpp
index fa772b69b..ed622c135 100644
--- a/src/geom/CircularArc.cpp
+++ b/src/geom/CircularArc.cpp
@@ -317,6 +317,16 @@ CircularArc::getDirectionPoint() const
return CircularArcs::getDirectionPoint(getCenter(), getRadius(), theta0(), getOrientation() == algorithm::Orientation::COUNTERCLOCKWISE);
}
+CoordinateXY
+CircularArc::getReverseDirectionPoint() const
+{
+ if (isLinear()) {
+ return p0();
+ }
+
+ return CircularArcs::getDirectionPoint(getCenter(), getRadius(), theta2(), getOrientation() != algorithm::Orientation::COUNTERCLOCKWISE);
+}
+
Envelope
CircularArc::getEnvelope() const
{
diff --git a/src/operation/overlayng/OverlayUtil.cpp b/src/operation/overlayng/OverlayUtil.cpp
index 077ca2efb..275ca422f 100644
--- a/src/operation/overlayng/OverlayUtil.cpp
+++ b/src/operation/overlayng/OverlayUtil.cpp
@@ -348,7 +348,7 @@ OverlayUtil::getDirectionPoint(const CoordinateSequence& pts, bool forward, bool
return arc.getDirectionPoint();
} else {
CircularArc arc(pts, pts.size() - 3);
- return algorithm::CircularArcs::getDirectionPoint(arc.getCenter(), arc.getRadius(), arc.theta2(), !arc.isCCW());
+ return arc.getReverseDirectionPoint();
}
}
commit fac5376f75356ad9cd8b6f08ef0336f1d56e006c
Author: Daniel Baston <dbaston at gmail.com>
Date: Mon Jun 8 08:35:59 2026 -0400
CircularArc: handle direction point of degenerate arcs
diff --git a/src/geom/CircularArc.cpp b/src/geom/CircularArc.cpp
index dd5a56d85..fa772b69b 100644
--- a/src/geom/CircularArc.cpp
+++ b/src/geom/CircularArc.cpp
@@ -310,6 +310,10 @@ CircularArc::getArea() const {
CoordinateXY
CircularArc::getDirectionPoint() const
{
+ if (isLinear()) {
+ return p2();
+ }
+
return CircularArcs::getDirectionPoint(getCenter(), getRadius(), theta0(), getOrientation() == algorithm::Orientation::COUNTERCLOCKWISE);
}
commit c06102e28158cd1393a68c767f36087b249e4ee0
Author: Daniel Baston <dbaston at gmail.com>
Date: Sun Jun 7 11:47:50 2026 -0400
CircularArcIntersector: Improve robustness of arc/segment endpoint intersection
diff --git a/include/geos/algorithm/CircularArcs.h b/include/geos/algorithm/CircularArcs.h
index a32a197cb..0aa868e23 100644
--- a/include/geos/algorithm/CircularArcs.h
+++ b/include/geos/algorithm/CircularArcs.h
@@ -65,6 +65,21 @@ public:
static void interpolateZM(const geom::CoordinateSequence &seq, size_t i0, const geom::CoordinateXY ¢er, bool isCCW, geom::CoordinateXY &pt, double
&z, double &m);
+ /** Determines whether and where a circle intersects a line.
+ *
+ * @param center The center point of the circle
+ * @param r The radius of the circle
+ * @param p0 One point defining a line of infinite length
+ * @param p1 Second point defining a line of infinite length
+ * @param isect0 Set to the first intersection point, if it exists
+ * @param isect1 Set to the second intersection point, if it exists
+ * @return The number of intersection points
+ */
+ static int
+ circleIntersectsLine(const geom::CoordinateXY& center, double r,
+ const geom::CoordinateXY& p0, const geom::CoordinateXY& p1,
+ geom::CoordinateXY& isect0, geom::CoordinateXY& isect1);
+
/** Determines whether and where a circle intersects a line segment.
*
* @param center The center point of the circle
diff --git a/src/algorithm/CircularArcIntersector.cpp b/src/algorithm/CircularArcIntersector.cpp
index c63cd2c55..ad79ee190 100644
--- a/src/algorithm/CircularArcIntersector.cpp
+++ b/src/algorithm/CircularArcIntersector.cpp
@@ -101,6 +101,23 @@ CircularArcIntersector::hasIntersection(const geom::CoordinateXY &p) const {
return false;
}
+static CoordinateXY&
+closestPoint(CoordinateXY& p0, CoordinateXY& p1, int n, const CoordinateXY& q)
+{
+ if (n < 2) {
+ return p0;
+ }
+
+ const double d0 = p0.distance(q);
+ const double d1 = p1.distance(q);
+
+ if (d0 < d1) {
+ return p0;
+ }
+
+ return p1;
+}
+
void
CircularArcIntersector::intersects(const CircularArc& arc, const CoordinateSequence& seq, std::size_t segPos0, std::size_t segPos1, bool useSegEndpoints)
{
@@ -117,14 +134,35 @@ CircularArcIntersector::intersects(const CircularArc& arc, const CoordinateSeque
const CoordinateXY& c = arc.getCenter();
const double r = arc.getRadius();
- CoordinateXYZM isect0, isect1;
- auto n = CircularArcs::circleIntersectsSegment(c, r, seq.getAt<CoordinateXY>(segPos0), seq.getAt<CoordinateXY>(segPos1), isect0, isect1);
+ CoordinateXY isect0, isect1;
+ const auto nPointsIntersectingLine = CircularArcs::circleIntersectsLine(c, r, seq.getAt<CoordinateXY>(segPos0), seq.getAt<CoordinateXY>(segPos1), isect0, isect1);
- if (n > 0 && arc.containsPointOnCircle(isect0)) {
+ if (nPointsIntersectingLine == 0) {
+ result = NO_INTERSECTION;
+ return;
+ }
+
+ // Check for exact endpoint-endpoint intersections
+ // If found, replace the computed intersection points with an exact endpoint
+ const CoordinateXY& ap0 = arc.p0<CoordinateXY>();
+ const CoordinateXY& ap2 = arc.p2<CoordinateXY>();
+ const CoordinateXY& bp0 = seq.getAt<CoordinateXY>(segPos0);
+ const CoordinateXY& bp1 = seq.getAt<CoordinateXY>(segPos1);
+
+ if (ap0 == bp0 || ap0 == bp1) {
+ closestPoint(isect0, isect1, nPointsIntersectingLine, ap0) = ap0;
+ }
+ if (ap2 == bp0 || ap2 == bp1) {
+ closestPoint(isect0, isect1, nPointsIntersectingLine, ap2) = ap2;
+ }
+
+ Envelope segEnv(bp0, bp1);
+
+ if (nPointsIntersectingLine > 0 && segEnv.contains(isect0) && arc.containsPointOnCircle(isect0)) {
addArcSegmentIntersectionPoint(isect0, arc, seq, segPos0, segPos1, useSegEndpoints);
}
- if (n > 1 && arc.containsPointOnCircle(isect1)) {
+ if (nPointsIntersectingLine > 1 && segEnv.contains(isect1) && arc.containsPointOnCircle(isect1)) {
addArcSegmentIntersectionPoint(isect1, arc, seq, segPos0, segPos1, useSegEndpoints);
}
diff --git a/src/algorithm/CircularArcs.cpp b/src/algorithm/CircularArcs.cpp
index eeb779e2b..07f2bd630 100644
--- a/src/algorithm/CircularArcs.cpp
+++ b/src/algorithm/CircularArcs.cpp
@@ -309,22 +309,26 @@ CircularArcs::expandEnvelope(geom::Envelope& e, const geom::CoordinateXY& p0, co
}
int
-CircularArcs::circleIntersectsSegment(const CoordinateXY& center, double r,
- const CoordinateXY& p0, const CoordinateXY& p1,
- CoordinateXY& ret0, CoordinateXY& ret1)
+CircularArcs::circleIntersectsLine(const geom::CoordinateXY& center, double r,
+ const geom::CoordinateXY& p0, const geom::CoordinateXY& p1,
+ geom::CoordinateXY& isect0, geom::CoordinateXY& isect1)
{
const double& x0 = center.x;
const double& y0 = center.y;
- Envelope segEnv(p0, p1);
-
- CoordinateXY isect0, isect1;
- int n = 0;
-
if (p1.x == p0.x) {
// vertical line
double x = p1.x;
+ if (p1.x > x0 + r || p1.x < x0 - r) {
+ return 0;
+ }
+
+ if (p1.x == x0 + r || p1.x == x0 - r) {
+ isect0 = {p1.x, y0};
+ return 1;
+ }
+
double A = 1;
double B = -2*y0;
double C = x*x - 2*x*x0 + x0*x0 + y0*y0 - r*r;
@@ -335,41 +339,73 @@ CircularArcs::circleIntersectsSegment(const CoordinateXY& center, double r,
isect0 = {x, Y1};
isect1 = {x, Y2};
+
+ return 2;
}
- else {
- double m = (p1.y - p0.y) / (p1.x - p0.x);
- double b = p1.y - p1.x*m;
- // Ax^2 + Bx + C = 0
- double A = 1 + m*m;
- double B = -2*x0 + 2*m*b - 2*m*y0;
- double C = x0*x0 + b*b - 2*b*y0 + y0*y0 - r*r;
+ const double m = (p1.y - p0.y) / (p1.x - p0.x);
+ const double b = p1.y - p1.x*m;
- double d = std::sqrt(B*B - 4*A*C);
- double X1 = (-B + d)/(2*A);
- double X2 = (-B - d)/(2*A);
+ // Take equation defining circle: (x - x0)^2 + (y - y0)^2 = r^2
+ // Substitute in equation defining line: y = mx + b
+ // Rearrange into standard quadratic form: Ax^2 + Bx + C = 0
+ const double A = 1 + m*m;
+ const double B = -2*x0 + 2*m*b - 2*m*y0;
+ const double C = x0*x0 + b*b - 2*b*y0 + y0*y0 - r*r;
- // TODO use robust quadratic solver such as https://github.com/archermarx/quadratic ?
- // auto [X1, X2] = quadratic::solve(A, B, C);
+ const double dd = B*B - 4*A*C;
- isect0 = {X1, m* X1 + b};
- isect1 = {X2, m* X2 + b};
+ if (dd < 0) {
+ return 0;
}
+ // TODO use robust quadratic solver such as https://github.com/archermarx/quadratic ?
+ // auto [X1, X2] = quadratic::solve(A, B, C);
+
+ const double d = std::sqrt(dd);
+ const double X1 = (-B + d)/(2*A);
+ const double X2 = (-B - d)/(2*A);
+
+ isect0 = {X1, m* X1 + b};
+ if (d == 0) {
+ return 1;
+ }
+
+ isect1 = {X2, m* X2 + b};
+ return 2;
+}
+
+int
+CircularArcs::circleIntersectsSegment(const CoordinateXY& center, double r,
+ const CoordinateXY& p0, const CoordinateXY& p1,
+ CoordinateXY& ret0, CoordinateXY& ret1)
+{
+ CoordinateXY isect0, isect1;
+
+ auto n = circleIntersectsLine(center, r, p0, p1, isect0, isect1);
+
+ if (n == 0) {
+ return 0;
+ }
+
+ Envelope segEnv(p0, p1);
+
if (segEnv.intersects(isect0)) {
ret0 = isect0;
- if (segEnv.intersects(isect1) && !isect1.equals2D(isect0)) {
+ if (n > 1 && segEnv.intersects(isect1)) {
ret1 = isect1;
- n = 2;
- } else {
- n = 1;
+ return 2;
}
- } else if (segEnv.intersects(isect1)) {
- ret0 = isect1;
- n = 1;
+
+ return 1;
}
- return n;
+ if (segEnv.intersects(isect1)) {
+ ret0 = isect1;
+ return 1;
+ }
+
+ return 0;
}
bool
diff --git a/tests/unit/algorithm/CircularArcIntersectorTest.cpp b/tests/unit/algorithm/CircularArcIntersectorTest.cpp
index 928a3df90..e9899213d 100644
--- a/tests/unit/algorithm/CircularArcIntersectorTest.cpp
+++ b/tests/unit/algorithm/CircularArcIntersectorTest.cpp
@@ -1712,6 +1712,26 @@ void object::test<85>()
ensure_equals(cai1.getPoint(0), cai2.getPoint(0));
}
+template<>
+template<>
+void object::test<86>()
+{
+ // arc/segment with single endpoint intersection
+
+ CircularArcIntersector cai;
+
+ auto arc = CircularArc::create(XY{-0.84333168960872817, -0.1241204385865409}, XY{-0.96461319189094762, -0.20729352391791989}, XY{-1.1087328289561522, -0.23655946199558511});
+
+ CoordinateSequence seg{
+ XY{-1.277492733105557, 0.28999497465537871}, XY{-0.84333168960872817, -0.1241204385865409}};
+
+ cai.intersects(arc, seg, 0, 1, false);
+
+ ensure_equals(cai.getNumPoints(), 1u);
+ // ensure endpoint intersection is represented exactly, not with distance() == 0
+ ensure_equals(cai.getPoint(0), seg.getAt<CoordinateXY>(1));
+}
+
// TODO: check Z values of arc result centerpoints
// TODO: add tests for seg/seg
commit f14c453ade3f86a98b3dafeaf5cca750593d4e70
Author: Daniel Baston <dbaston at gmail.com>
Date: Fri May 15 14:12:27 2026 -0400
Edge: simplify direction() logic
diff --git a/include/geos/operation/overlayng/Edge.h b/include/geos/operation/overlayng/Edge.h
index 83f67dc6a..b88a44fe7 100644
--- a/include/geos/operation/overlayng/Edge.h
+++ b/include/geos/operation/overlayng/Edge.h
@@ -230,25 +230,22 @@ public:
throw util::GEOSException("Edge must have >= 2 points");
}
- const geom::Coordinate& p0 = pts->getAt(0);
- const geom::Coordinate& p1 = pts->getAt(1);
- const geom::Coordinate& pn0 = pts->getAt(pts->size() - 1);
- const geom::Coordinate& pn1 = pts->getAt(pts->size() - 2);
+ const geom::CoordinateXY& p0 = getCoordinate(0);
+ const geom::CoordinateXY& pn0 = getCoordinate(pts->size() - 1);
- int cmp = 0;
- int cmp0 = p0.compareTo(pn0);
- if (cmp0 != 0) cmp = cmp0;
-
- if (cmp == 0) {
- int cmp1 = p1.compareTo(pn1);
- if (cmp1 != 0) cmp = cmp1;
+ const int cmp0 = p0.compareTo(pn0);
+ if (cmp0 != 0) {
+ return cmp0 == -1;
}
- if (cmp == 0) {
- throw util::GEOSException("Edge direction cannot be determined because endpoints are equal");
+ const geom::CoordinateXY& p1 = getCoordinate(1);
+ const geom::CoordinateXY& pn1 = getCoordinate(pts->size() - 2);
+ const int cmp1 = p1.compareTo(pn1);
+ if (cmp1 != 0) {
+ return cmp1 == -1;
}
- return cmp == -1;
+ throw util::GEOSException("Edge direction cannot be determined because endpoints are equal");
};
/**
commit dc3935ba35f5fd6328bddb3637c7e0eeba307dd8
Author: Daniel Baston <dbaston at gmail.com>
Date: Fri May 15 14:09:56 2026 -0400
NodableArcString: Split ArcStrings when intersection point is already a node
diff --git a/src/noding/NodableArcString.cpp b/src/noding/NodableArcString.cpp
index d16acafd5..4b824b3fa 100644
--- a/src/noding/NodableArcString.cpp
+++ b/src/noding/NodableArcString.cpp
@@ -16,6 +16,7 @@
#include <geos/algorithm/Angle.h>
#include <algorithm>
+#include <iomanip>
using geos::geom::CoordinateXYZM;
using geos::geom::CircularArc;
@@ -26,7 +27,7 @@ static double
pseudoAngleDiffCCW(double paStart, double pa) {
double diff = pa - paStart;
- if (diff <= 0) {
+ if (diff < 0) {
diff += 4;
}
@@ -45,31 +46,32 @@ std::vector<CoordinateXYZM> prepareArcPoints(const CircularArc& arc, std::vector
{
const bool isCCW = arc.getOrientation() == algorithm::Orientation::COUNTERCLOCKWISE;
const geom::CoordinateXY& center = arc.getCenter();
- const double paStart = geom::Quadrant::pseudoAngle(center, arc.p0());
+
+ // Some potential split points may be skipped, for example, if they would create an arc section that is
+ // too short to have a constructed midpoint. Because the results of this logic could depend on the
+ // direction in which the arc is processed, we reverse clockwise arcs and then reverse the list of
+ // retained points.
+ const double paStart = geom::Quadrant::pseudoAngle(center, isCCW ? arc.p0() : arc.p2());
std::vector<CoordinateXYZM> retained;
// Add starting point of input arc
{
CoordinateXYZM p0;
- arc.getCoordinateSequence()->getAt(arc.getCoordinatePosition(), p0);
+ arc.getCoordinateSequence()->getAt(arc.getCoordinatePosition() + (isCCW ? 0 : 2), p0);
retained.push_back(p0);
}
- std::sort(splitPoints.begin(), splitPoints.end(), [¢er, paStart, isCCW](const auto& p0, const auto& p1) {
+ std::sort(splitPoints.begin(), splitPoints.end(), [¢er, paStart](const auto& p0, const auto& p1) {
double pa0 = geom::Quadrant::pseudoAngle(center, p0);
double pa1 = geom::Quadrant::pseudoAngle(center, p1);
- if (isCCW) {
- return pseudoAngleDiffCCW(paStart, pa0) < pseudoAngleDiffCCW(paStart, pa1);
- } else {
- return pseudoAngleDiffCCW(paStart, pa0) > pseudoAngleDiffCCW(paStart, pa1);
- }
+ return pseudoAngleDiffCCW(paStart, pa0) < pseudoAngleDiffCCW(paStart, pa1);
});
// Add ending point of input arc
{
CoordinateXYZM p2;
- arc.getCoordinateSequence()->getAt(arc.getCoordinatePosition() + 2, p2);
+ arc.getCoordinateSequence()->getAt(arc.getCoordinatePosition() + (isCCW ? 2 : 0), p2);
splitPoints.push_back(p2);
}
@@ -90,23 +92,21 @@ std::vector<CoordinateXYZM> prepareArcPoints(const CircularArc& arc, std::vector
continue;
}
- const auto p1 = algorithm::CircularArcs::getMidpoint(p0, p2, center, arc.getRadius(), isCCW);
+ // Calculate the midpoint of an arc between p0 and p2.
+ // We don't actually use the calculated point here, we just want to make sure that
+ // the arc from p0 to p2 is long enough to contain a midpoint.
+ const geom::CoordinateXY p1 = algorithm::CircularArcs::getMidpoint(p0, p2, center, arc.getRadius(), true);
if (p1.equals2D(p0) || p1.equals2D(p2)) {
continue;
}
- // Reject split point where sub-arc midpoint doesn't fall inside arc
- if (!arc.containsPointOnCircle(p1)) {
- continue;
- }
-
// Reject split point where computed doesn't fall between endpoints
- const double t0 = algorithm::Angle::normalizePositive(arc.theta0());
+ const double t0 = algorithm::Angle::normalizePositive(isCCW ? arc.theta0() : arc.theta2());
const double t1 = algorithm::Angle::normalizePositive(algorithm::CircularArcs::getAngle(p1, center));
- const double t2 = algorithm::Angle::normalizePositive(algorithm::CircularArcs::getAngle(p2, center));
+ const double t2 = algorithm::Angle::normalizePositive(isCCW ? arc.theta2() : arc.theta0());
- if (algorithm::Angle::isWithinCCW(t1, t0, t2) != isCCW) {
+ if (!algorithm::Angle::isWithinCCW(t1, t0, t2)) {
continue;
}
@@ -116,7 +116,7 @@ std::vector<CoordinateXYZM> prepareArcPoints(const CircularArc& arc, std::vector
// Make sure that endpoint of split arc is unchanged
{
CoordinateXYZM p2;
- arc.getCoordinateSequence()->getAt(arc.getCoordinatePosition() + 2, p2);
+ arc.getCoordinateSequence()->getAt(arc.getCoordinatePosition() + (isCCW ? 2 : 0), p2);
CoordinateXYZM& back = retained.back();
if (!back.equals2D(p2)) {
@@ -131,6 +131,10 @@ std::vector<CoordinateXYZM> prepareArcPoints(const CircularArc& arc, std::vector
}
}
+ if (!isCCW) {
+ std::reverse(retained.begin(), retained.end());
+ }
+
return retained;
}
@@ -146,16 +150,19 @@ NodableArcString::getNoded(std::vector<std::unique_ptr<ArcString>>& splitArcs) {
const int orientation = toSplit.getOrientation();
bool arcIsSplit = true;
+ bool createArcString = true;
const bool preserveControlPoint = true;
std::vector<CoordinateXYZM> arcPoints;
const auto it = m_adds.find(arcIndex);
if (it == m_adds.end()) {
arcIsSplit = false;
+ createArcString = false;
} else {
arcPoints = prepareArcPoints(toSplit, it->second);
if (arcPoints.size() == 2) {
// All added nodes collapsed
+ // Still need to know if arc was split.
arcIsSplit = false;
}
}
@@ -168,6 +175,12 @@ NodableArcString::getNoded(std::vector<std::unique_ptr<ArcString>>& splitArcs) {
std::size_t dstPos = dstSeq->getSize() - 3;
arcs.emplace_back(*dstSeq, dstPos, center, radius, orientation);
+ if (createArcString) {
+ splitArcs.push_back(std::make_unique<NodableArcString>(std::move(arcs), std::move(dstSeq), m_constructZ, m_constructM, getData()));
+ dstSeq = std::make_unique<geom::CoordinateSequence>(0, m_constructZ, m_constructM);
+ arcs.clear();
+ }
+
continue;
}
diff --git a/tests/unit/capi/GEOSNodeTest.cpp b/tests/unit/capi/GEOSNodeTest.cpp
index 8feb5df96..63d893a55 100644
--- a/tests/unit/capi/GEOSNodeTest.cpp
+++ b/tests/unit/capi/GEOSNodeTest.cpp
@@ -235,10 +235,10 @@ void object::test<10>()
ensure(result_ != nullptr);
expected_ = fromWKT("MULTICURVE ZM ("
- "CIRCULARSTRING ZM (-1 0 3 4, -0.7071067812 0.7071067812 2.5 6.5, -5.5511151231e-17 1 2 9),"
- "CIRCULARSTRING ZM (-5.5511151231e-17 1 2 9, 0.7071067812 0.7071067812 3 8, 1 0 4 7),"
- "CIRCULARSTRING ZM (-1 2 NaN 9, -0.2928932188 1.7071067812 NaN 9, -5.5511151231e-17 1 2 9),"
- "CIRCULARSTRING ZM (-5.5511151231e-17 1 2 9, -0.2928932188 0.2928932188 2.5 6.5, -1 0 3 4))");
+ "CIRCULARSTRING ZM (-1 0 3 4, -0.7071067812 0.7071067812 2.5 6.5, 0 1 2 9),"
+ "CIRCULARSTRING ZM (0 1 2 9, 0.7071067812 0.7071067812 3 8, 1 0 4 7),"
+ "CIRCULARSTRING ZM (-1 2 NaN 9, -0.2928932188 1.7071067812 NaN 9, 0 1 2 9),"
+ "CIRCULARSTRING ZM (0 1 2 9, -0.2928932188 0.2928932188 2.5 13, -1 0 3 17))");
ensure_equals(GEOSGetNumGeometries(result_), 4);
diff --git a/tests/unit/geom/CircularStringTest.cpp b/tests/unit/geom/CircularStringTest.cpp
index ee4e26058..3678746c9 100644
--- a/tests/unit/geom/CircularStringTest.cpp
+++ b/tests/unit/geom/CircularStringTest.cpp
@@ -138,9 +138,9 @@ void object::test<3>()
// Overlay
ensure_THROW(cs_->Union(), geos::util::UnsupportedOperationException);
- ensure_equals_geometry(cs_->Union(cs_.get()).get(), static_cast<const Geometry*>(cs_.get()));
+ ensure_equals(cs_->Union(cs_.get())->getLength(), cs_->getLength());
ensure(cs_->difference(cs_.get())->isEmpty());
- ensure_equals_geometry(cs_->intersection(cs_.get()).get(), static_cast<const Geometry*>(cs_.get()));
+ ensure_equals(cs_->intersection(cs_.get())->getLength(), cs_->getLength());
ensure(cs_->symDifference(cs_.get())->isEmpty());
// Distance
diff --git a/tests/unit/geom/CurvePolygonTest.cpp b/tests/unit/geom/CurvePolygonTest.cpp
index 07c3faad1..012fa8305 100644
--- a/tests/unit/geom/CurvePolygonTest.cpp
+++ b/tests/unit/geom/CurvePolygonTest.cpp
@@ -166,9 +166,9 @@ void object::test<3>()
// Overlay
ensure_THROW(cp_->Union(), geos::util::UnsupportedOperationException);
- ensure_equals_geometry(cp_->Union(cp_.get()).get(), static_cast<const Geometry*>(cp_.get()));
+ ensure_equals(cp_->Union(cp_.get())->getLength(), cp_->getLength());
ensure(cp_->difference(cp_.get())->isEmpty());
- ensure_equals_geometry(cp_->intersection(cp_.get()).get(), static_cast<const Geometry*>(cp_.get()));
+ ensure_equals(cp_->intersection(cp_.get())->getLength(), cp_->getLength());
ensure(cp_->symDifference(cp_.get())->isEmpty());
// Distance
diff --git a/tests/unit/operation/overlayng/OverlayNGTest.cpp b/tests/unit/operation/overlayng/OverlayNGTest.cpp
index d7077e57a..46b8155fa 100644
--- a/tests/unit/operation/overlayng/OverlayNGTest.cpp
+++ b/tests/unit/operation/overlayng/OverlayNGTest.cpp
@@ -780,7 +780,7 @@ template<>
template<>
void object::test<60>()
{
- set_test_name("MultiSurface / MultiPoint -> MultiSurface");
+ set_test_name("Union of MultiSurface / MultiPoint -> MultiSurface");
std::string a = "MULTISURFACE(CURVEPOLYGON (COMPOUNDCURVE( CIRCULARSTRING(-5 0, 0 5, 5 0), (5 0, -5 0))), ((4 4, 5 4, 5 5, 4 4)))";
std::string b = "MULTIPOINT (0 0, 4 3, 2 2, 4.5 4.1)";
@@ -789,4 +789,33 @@ void object::test<60>()
testOverlay(a, b, exp, OverlayNG::UNION, 0);
}
+template<>
+template<>
+void object::test<61>()
+{
+ set_test_name("Union of CurvePolygon / CurvePolygon -> CurvePolygon");
+
+ std::string a = "CURVEPOLYGON (CIRCULARSTRING (-5 0, 0 5, 5 0, 0 4, -5 0))";
+ std::string b = "CURVEPOLYGON (COMPOUNDCURVE((-5 0, 5 0), CIRCULARSTRING (5 0, 0 4, -5 0)))";
+
+ std::string exp = "CURVEPOLYGON (COMPOUNDCURVE (CIRCULARSTRING (-5 0, 0 5, 5 0), (5 0, -5 0)))";
+
+ testOverlay(a, b, exp, OverlayNG::UNION, 0);
+}
+
+template<>
+template<>
+void object::test<62>()
+{
+ set_test_name("CircularString self-union");
+
+ std::string a = "CIRCULARSTRING (0 0, 1 1, 2 0, 3 -1, 4 0)";
+
+ // Noding causes the CircularString to split at (2, 0)
+ // OverlayNG does not merge output lines, so we get a MultiCurve.
+ std::string exp = "MULTICURVE (CIRCULARSTRING (0 0, 1 1, 2 0), CIRCULARSTRING (2 0, 3 -1, 4 0))";
+
+ testOverlay(a, a, exp, OverlayNG::UNION, 0);
+}
+
} // namespace tut
commit 1aeef7493bc6d0804d094021fc1013a7227791af
Author: Daniel Baston <dbaston at gmail.com>
Date: Thu May 14 15:57:06 2026 -0400
CircularArcIntersector: Calc intersection points independently of argument order
diff --git a/src/algorithm/CircularArcIntersector.cpp b/src/algorithm/CircularArcIntersector.cpp
index b5e6adcc3..c63cd2c55 100644
--- a/src/algorithm/CircularArcIntersector.cpp
+++ b/src/algorithm/CircularArcIntersector.cpp
@@ -163,11 +163,15 @@ CircularArcIntersector::intersects(const CircularArc& arc1, const CircularArc& a
reset();
- const auto& c1 = arc1.getCenter();
- const auto& c2 = arc2.getCenter();
+ // Normalize arguments such that the computed intersection points do not depend
+ // on the order of the input arcs
+ const bool swapArgs = arc1.getCenter().compareTo(arc2.getCenter()) > 0;
- const auto r1 = arc1.getRadius();
- const auto r2 = arc2.getRadius();
+ const auto c1 = swapArgs ? arc2.getCenter() : arc1.getCenter();
+ const auto c2 = swapArgs ? arc1.getCenter() : arc2.getCenter();
+
+ const auto r1 = swapArgs ? arc2.getRadius() : arc1.getRadius();
+ const auto r2 = swapArgs ? arc1.getRadius() : arc2.getRadius();
const auto d = c1.distance(c2);
diff --git a/tests/unit/algorithm/CircularArcIntersectorTest.cpp b/tests/unit/algorithm/CircularArcIntersectorTest.cpp
index 6ad649d02..928a3df90 100644
--- a/tests/unit/algorithm/CircularArcIntersectorTest.cpp
+++ b/tests/unit/algorithm/CircularArcIntersectorTest.cpp
@@ -1692,6 +1692,26 @@ void object::test<84>()
ensure_equals(cai.getPoint(0), XY{31.8, 68.2});
}
+template<>
+template<>
+void object::test<85>()
+{
+ set_test_name("computed test points are independent of argument order");
+
+ auto arc0 = CircularArc::create(XY{0, 5}, XY{3, 4}, XY{4, 3});
+ auto arc1 = CircularArc::create(XY{3, 5}, XY{4, 4}, XY{2, 3});
+
+ CircularArcIntersector cai1;
+ cai1.intersects(arc0, arc1);
+ ensure_equals(cai1.getNumPoints(), 1u);
+
+ CircularArcIntersector cai2;
+ cai2.intersects(arc1, arc0);
+ ensure_equals(cai2.getNumPoints(), 1u);
+
+ ensure_equals(cai1.getPoint(0), cai2.getPoint(0));
+}
+
// TODO: check Z values of arc result centerpoints
// TODO: add tests for seg/seg
commit b0506ea14bdc0339fb8ab61eb95ba5ec26762332
Author: Daniel Baston <dbaston at gmail.com>
Date: Thu May 14 13:51:20 2026 -0400
OverlayUtil: Add getDirectionPoint
diff --git a/include/geos/operation/overlayng/OverlayUtil.h b/include/geos/operation/overlayng/OverlayUtil.h
index 25fa99312..5a83bda6a 100644
--- a/include/geos/operation/overlayng/OverlayUtil.h
+++ b/include/geos/operation/overlayng/OverlayUtil.h
@@ -205,6 +205,9 @@ public:
*/
static bool round(const Point* pt, const PrecisionModel* pm, Coordinate& rsltCoord);
+ static geom::CoordinateXY
+ getDirectionPoint(const CoordinateSequence& pts, bool forward, bool isCurved);
+
template<typename T>
static void moveGeometry(std::vector<std::unique_ptr<T>>& inGeoms, std::vector<std::unique_ptr<Geometry>>& outGeoms)
{
diff --git a/src/operation/overlayng/OverlayGraph.cpp b/src/operation/overlayng/OverlayGraph.cpp
index 7fd39fefd..8fc1de1ca 100644
--- a/src/operation/overlayng/OverlayGraph.cpp
+++ b/src/operation/overlayng/OverlayGraph.cpp
@@ -14,8 +14,8 @@
#include <geos/operation/overlayng/OverlayGraph.h>
-#include <geos/algorithm/CircularArcs.h>
#include <geos/operation/overlayng/Edge.h>
+#include <geos/operation/overlayng/OverlayUtil.h>
#include <geos/geom/Coordinate.h>
#include <geos/geom/CoordinateSequence.h>
#include <geos/geom/CircularArc.h>
@@ -105,28 +105,6 @@ OverlayGraph::createEdgePair(const std::shared_ptr<const CoordinateSequence>& pt
return e0;
}
-static CoordinateXY
-getDirectionPoint(const CoordinateSequence& pts, bool forward, bool isCurved)
-{
- if (isCurved) {
- assert(pts.size() >= 3);
- if (forward) {
- CircularArc arc(pts, 0);
- return arc.getDirectionPoint();
- } else {
- CircularArc arc(pts, pts.size() - 3);
- return algorithm::CircularArcs::getDirectionPoint(arc.getCenter(), arc.getRadius(), arc.theta2(), !arc.isCCW());
- }
- }
-
- assert(pts.size() >= 2);
- if (forward) {
- return pts.getAt<CoordinateXY>(1);
- }
-
- return pts.getAt<CoordinateXY>(pts.size() - 2);
-}
-
/*private*/
OverlayEdge*
OverlayGraph::createOverlayEdge(const std::shared_ptr<const CoordinateSequence>& pts, OverlayLabel* lbl, bool direction, bool isCurved)
@@ -134,7 +112,7 @@ OverlayGraph::createOverlayEdge(const std::shared_ptr<const CoordinateSequence>&
assert(pts->size() >= 2);
CoordinateXYZM origin;
- const CoordinateXY dirPt = getDirectionPoint(*pts, direction, isCurved);
+ const CoordinateXY dirPt = OverlayUtil::getDirectionPoint(*pts, direction, isCurved);
if (direction) {
pts->getAt(0, origin);
diff --git a/src/operation/overlayng/OverlayUtil.cpp b/src/operation/overlayng/OverlayUtil.cpp
index 313837a34..077ca2efb 100644
--- a/src/operation/overlayng/OverlayUtil.cpp
+++ b/src/operation/overlayng/OverlayUtil.cpp
@@ -338,6 +338,28 @@ OverlayUtil::round(const Point* pt, const PrecisionModel* pm, Coordinate& rsltCo
return true;
}
+CoordinateXY
+OverlayUtil::getDirectionPoint(const CoordinateSequence& pts, bool forward, bool isCurved)
+{
+ if (isCurved) {
+ assert(pts.size() >= 3);
+ if (forward) {
+ CircularArc arc(pts, 0);
+ return arc.getDirectionPoint();
+ } else {
+ CircularArc arc(pts, pts.size() - 3);
+ return algorithm::CircularArcs::getDirectionPoint(arc.getCenter(), arc.getRadius(), arc.theta2(), !arc.isCCW());
+ }
+ }
+
+ assert(pts.size() >= 2);
+ if (forward) {
+ return pts.getAt<CoordinateXY>(1);
+ }
+
+ return pts.getAt<CoordinateXY>(pts.size() - 2);
+}
+
} // namespace geos.operation.overlayng
commit 43edb6d35a6c66fe82712033bb3abd53f89fce3a
Author: Daniel Baston <dbaston at gmail.com>
Date: Tue Apr 14 06:12:17 2026 -0400
EdgeNodingBuilder: Automatically use SimpleNoder for curved inputs
diff --git a/include/geos/operation/overlayng/EdgeNodingBuilder.h b/include/geos/operation/overlayng/EdgeNodingBuilder.h
index d981a51d5..8d34ceffd 100644
--- a/include/geos/operation/overlayng/EdgeNodingBuilder.h
+++ b/include/geos/operation/overlayng/EdgeNodingBuilder.h
@@ -18,6 +18,7 @@
#pragma once
+#include <geos/algorithm/CircularArcIntersector.h>
#include <geos/algorithm/LineIntersector.h>
#include <geos/geom/CoordinateSequence.h>
#include <geos/geom/CircularArc.h>
@@ -28,6 +29,7 @@
#include <geos/geom/LineString.h>
#include <geos/geom/Polygon.h>
#include <geos/geom/SimpleCurve.h>
+#include <geos/noding/ArcIntersectionAdder.h>
#include <geos/noding/IntersectionAdder.h>
#include <geos/noding/Noder.h>
#include <geos/noding/SegmentString.h>
@@ -95,7 +97,9 @@ private:
// For use in createFloatingPrecisionNoder()
algorithm::LineIntersector lineInt;
+ algorithm::CircularArcIntersector arcInt;
noding::IntersectionAdder intAdder;
+ noding::ArcIntersectionAdder arcIntAdder;
std::unique_ptr<noding::Noder> internalNoder;
std::unique_ptr<noding::Noder> spareInternalNoder;
// EdgeSourceInfo*, Edge* owned by EdgeNodingBuilder, stored in deque
@@ -103,6 +107,7 @@ private:
std::deque<Edge> edgeQue;
bool inputHasZ;
bool inputHasM;
+ bool inputHasCurves;
/**
* Gets a noder appropriate for the precision model supplied.
@@ -207,8 +212,10 @@ public:
, hasEdges{{false,false}}
, clipEnv(nullptr)
, intAdder(lineInt)
+ , arcIntAdder(arcInt)
, inputHasZ(false)
, inputHasM(false)
+ , inputHasCurves(false)
{};
~EdgeNodingBuilder();
diff --git a/src/operation/overlayng/EdgeNodingBuilder.cpp b/src/operation/overlayng/EdgeNodingBuilder.cpp
index 13c123315..ded0c128b 100644
--- a/src/operation/overlayng/EdgeNodingBuilder.cpp
+++ b/src/operation/overlayng/EdgeNodingBuilder.cpp
@@ -23,6 +23,7 @@
#include <geos/noding/NodedSegmentString.h>
#include <geos/noding/MCIndexNoder.h>
#include <geos/noding/NodableArcString.h>
+#include <geos/noding/SimpleNoder.h>
#include <geos/noding/snapround/SnapRoundingNoder.h>
#include <geos/operation/overlayng/EdgeNodingBuilder.h>
#include <geos/operation/overlayng/EdgeMerger.h>
@@ -43,6 +44,7 @@ using geos::noding::NodableArcString;
using geos::noding::MCIndexNoder;
using geos::noding::ValidatingNoder;
using geos::noding::SegmentString;
+using geos::noding::SimpleNoder;
using geos::noding::NodedSegmentString;
using geos::noding::PathString;
@@ -83,16 +85,26 @@ EdgeNodingBuilder::createFixedPrecisionNoder(const PrecisionModel* p_pm)
std::unique_ptr<Noder>
EdgeNodingBuilder::createFloatingPrecisionNoder(bool doValidation)
{
- std::unique_ptr<MCIndexNoder> mcNoder(new MCIndexNoder());
- mcNoder->setSegmentIntersector(&intAdder);
+ std::unique_ptr<Noder> ret;
- if (doValidation) {
- spareInternalNoder = std::move(mcNoder);
+ if (inputHasCurves) {
+ auto simpleNoder = std::make_unique<SimpleNoder>();
+ simpleNoder->setArcIntersector(arcIntAdder);
+ ret = std::move(simpleNoder);
+ } else {
+ auto mcNoder = std::make_unique<MCIndexNoder>();
+ mcNoder->setSegmentIntersector(&intAdder);
+ ret = std::move(mcNoder);
+ }
+
+ // TODO: Support arcs in ValidatingNoder
+ if (doValidation && !inputHasCurves) {
+ spareInternalNoder = std::move(ret);
std::unique_ptr<Noder> validNoder(new ValidatingNoder(*spareInternalNoder));
return validNoder;
}
- return std::unique_ptr<Noder>(mcNoder.release());
+ return ret;
}
/*public*/
@@ -110,6 +122,7 @@ EdgeNodingBuilder::build(const Geometry* geom0, const Geometry* geom1)
{
inputHasZ = geom0->hasZ() || (geom1 != nullptr && geom1->hasZ());
inputHasM = geom0->hasM() || (geom1 != nullptr && geom1->hasM());
+ inputHasCurves = geom0->hasCurvedComponents() || (geom1 != nullptr && geom1->hasCurvedComponents());
add(geom0, 0);
add(geom1, 1);
diff --git a/src/operation/overlayng/OverlayNGRobust.cpp b/src/operation/overlayng/OverlayNGRobust.cpp
index e9af20896..258a7ecda 100644
--- a/src/operation/overlayng/OverlayNGRobust.cpp
+++ b/src/operation/overlayng/OverlayNGRobust.cpp
@@ -83,9 +83,6 @@ OverlayNGRobust::Union(const Geometry* a)
std::unique_ptr<Geometry>
OverlayNGRobust::Overlay(const Geometry* geom0, const Geometry* geom1, int opCode)
{
- geos::util::ensureNoCurvedComponents(geom0);
- geos::util::ensureNoCurvedComponents(geom1);
-
std::unique_ptr<Geometry> result;
std::runtime_error exOriginal("");
diff --git a/tests/unit/capi/GEOSDifferenceTest.cpp b/tests/unit/capi/GEOSDifferenceTest.cpp
index d739b0936..841f49188 100644
--- a/tests/unit/capi/GEOSDifferenceTest.cpp
+++ b/tests/unit/capi/GEOSDifferenceTest.cpp
@@ -63,6 +63,8 @@ template<>
template<>
void object::test<3>()
{
+ set_test_name("curved inputs");
+
geom1_ = fromWKT("CIRCULARSTRING (0 0, 1 1, 2 0)");
geom2_ = fromWKT("LINESTRING (1 0, 2 1)");
@@ -70,7 +72,11 @@ void object::test<3>()
ensure(geom2_);
result_ = GEOSDifference(geom1_, geom2_);
- ensure("curved geometry not supported", result_ == nullptr);
+ ensure(result_);
+
+ expected_ = fromWKT("MULTICURVE (CIRCULARSTRING (1.7071067812 0.7071067812, 1.9238795325 0.3826834324, 2 0), CIRCULARSTRING (0 0, 0.6173165676 0.9238795325, 1.7071067812 0.7071067812))");
+
+ ensure_geometry_equals(result_, expected_, 1e-8);
}
} // namespace tut
diff --git a/tests/unit/capi/GEOSIntersectionTest.cpp b/tests/unit/capi/GEOSIntersectionTest.cpp
index 4a03845df..a4aabb439 100644
--- a/tests/unit/capi/GEOSIntersectionTest.cpp
+++ b/tests/unit/capi/GEOSIntersectionTest.cpp
@@ -149,14 +149,16 @@ template<>
template<>
void object::test<8>()
{
+ set_test_name("curved input");
+
geom1_ = fromWKT("CIRCULARSTRING (0 0, 1 1, 2 0)");
geom2_ = fromWKT("LINESTRING (1 0, 2 1)");
- ensure(geom1_);
- ensure(geom2_);
-
result_ = GEOSIntersection(geom1_, geom2_);
- ensure("curved geometry not supported", result_ == nullptr);
+ ensure(result_);
+
+ expected_ = fromWKT("POINT (1.7071067812 0.7071067812)");
+ ensure_geometry_equals(result_, expected_, 1e-8);
}
// https://github.com/libgeos/geos/issues/1074
diff --git a/tests/unit/capi/GEOSSymDifferenceTest.cpp b/tests/unit/capi/GEOSSymDifferenceTest.cpp
index def092b68..63436e135 100644
--- a/tests/unit/capi/GEOSSymDifferenceTest.cpp
+++ b/tests/unit/capi/GEOSSymDifferenceTest.cpp
@@ -33,14 +33,16 @@ template<>
template<>
void object::test<2>()
{
+ set_test_name("curved inputs");
+
geom1_ = fromWKT("CIRCULARSTRING (0 0, 1 1, 2 0)");
geom2_ = fromWKT("LINESTRING (1 0, 2 1)");
- ensure(geom1_);
- ensure(geom2_);
-
result_ = GEOSSymDifference(geom1_, geom2_);
- ensure("curved geometry not supported", result_ == nullptr);
+ ensure(result_);
+
+ expected_ = fromWKT("MULTICURVE (CIRCULARSTRING (1.7071067812 0.7071067812, 1.9238795325 0.3826834324, 2 0), CIRCULARSTRING (0 0, 0.6173165676 0.9238795325, 1.7071067812 0.7071067812), (1.7071067812 0.7071067812, 2 1), (1 0, 1.7071067812 0.7071067812))");
+ ensure_geometry_equals(result_, expected_, 1e-8);
}
} // namespace tut
diff --git a/tests/unit/capi/GEOSUnionTest.cpp b/tests/unit/capi/GEOSUnionTest.cpp
index 68c90867b..4c7e89e29 100644
--- a/tests/unit/capi/GEOSUnionTest.cpp
+++ b/tests/unit/capi/GEOSUnionTest.cpp
@@ -66,14 +66,16 @@ template<>
template<>
void object::test<3>()
{
+ set_test_name("curved inputs");
+
geom1_ = fromWKT("CIRCULARSTRING (0 0, 1 1, 2 0)");
geom2_ = fromWKT("LINESTRING (1 0, 2 1)");
- ensure(geom1_);
- ensure(geom2_);
-
result_ = GEOSUnion(geom1_, geom2_);
- ensure("curved geometry not supported", result_ == nullptr);
+ ensure(result_);
+
+ expected_ = fromWKT("MULTICURVE (CIRCULARSTRING (1.7071067812 0.7071067812, 1.9238795325 0.3826834324, 2 0), CIRCULARSTRING (0 0, 0.6173165676 0.9238795325, 1.7071067812 0.7071067812), (1.7071067812 0.7071067812, 2 1), (1 0, 1.7071067812 0.7071067812))");
+ ensure_geometry_equals(result_, expected_, 1e-8);
}
} // namespace tut
diff --git a/tests/unit/geom/CircularStringTest.cpp b/tests/unit/geom/CircularStringTest.cpp
index a31b70769..ee4e26058 100644
--- a/tests/unit/geom/CircularStringTest.cpp
+++ b/tests/unit/geom/CircularStringTest.cpp
@@ -138,10 +138,10 @@ void object::test<3>()
// Overlay
ensure_THROW(cs_->Union(), geos::util::UnsupportedOperationException);
- ensure_THROW(cs_->Union(cs_.get()), geos::util::UnsupportedOperationException);
- ensure_THROW(cs_->difference(cs_.get()), geos::util::UnsupportedOperationException);
- ensure_THROW(cs_->intersection(cs_.get()), geos::util::UnsupportedOperationException);
- ensure_THROW(cs_->symDifference(cs_.get()), geos::util::UnsupportedOperationException);
+ ensure_equals_geometry(cs_->Union(cs_.get()).get(), static_cast<const Geometry*>(cs_.get()));
+ ensure(cs_->difference(cs_.get())->isEmpty());
+ ensure_equals_geometry(cs_->intersection(cs_.get()).get(), static_cast<const Geometry*>(cs_.get()));
+ ensure(cs_->symDifference(cs_.get())->isEmpty());
// Distance
ensure_equals(cs_->distance(cs_.get()), 0);
diff --git a/tests/unit/geom/CompoundCurveTest.cpp b/tests/unit/geom/CompoundCurveTest.cpp
index cc2f79cd8..4d52c66ef 100644
--- a/tests/unit/geom/CompoundCurveTest.cpp
+++ b/tests/unit/geom/CompoundCurveTest.cpp
@@ -164,10 +164,11 @@ void object::test<3>()
// Overlay
ensure_THROW(cc_->Union(), geos::util::UnsupportedOperationException);
- ensure_THROW(cc_->Union(cc_.get()), geos::util::UnsupportedOperationException);
- ensure_THROW(cc_->difference(cc_.get()), geos::util::UnsupportedOperationException);
- ensure_THROW(cc_->intersection(cc_.get()), geos::util::UnsupportedOperationException);
- ensure_THROW(cc_->symDifference(cc_.get()), geos::util::UnsupportedOperationException);
+ // TODO: Prevent overlay from degrading CompoundCurve into MultiCurve
+ //ensure_equals_geometry(cc_->Union(cc_.get()).get(), static_cast<const Geometry*>(cc_.get()));
+ ensure(cc_->difference(cc_.get())->isEmpty());
+ //ensure_equals_geometry(cc_->intersection(cc_.get()).get(), static_cast<const Geometry*>(cc_.get()));
+ ensure(cc_->symDifference(cc_.get())->isEmpty());
// Distance
ensure_equals(cc_->distance(cc_.get()), 0);
diff --git a/tests/unit/geom/CurvePolygonTest.cpp b/tests/unit/geom/CurvePolygonTest.cpp
index c41dab112..07c3faad1 100644
--- a/tests/unit/geom/CurvePolygonTest.cpp
+++ b/tests/unit/geom/CurvePolygonTest.cpp
@@ -166,10 +166,10 @@ void object::test<3>()
// Overlay
ensure_THROW(cp_->Union(), geos::util::UnsupportedOperationException);
- ensure_THROW(cp_->Union(cp_.get()), geos::util::UnsupportedOperationException);
- ensure_THROW(cp_->difference(cp_.get()), geos::util::UnsupportedOperationException);
- ensure_THROW(cp_->intersection(cp_.get()), geos::util::UnsupportedOperationException);
- ensure_THROW(cp_->symDifference(cp_.get()), geos::util::UnsupportedOperationException);
+ ensure_equals_geometry(cp_->Union(cp_.get()).get(), static_cast<const Geometry*>(cp_.get()));
+ ensure(cp_->difference(cp_.get())->isEmpty());
+ ensure_equals_geometry(cp_->intersection(cp_.get()).get(), static_cast<const Geometry*>(cp_.get()));
+ ensure(cp_->symDifference(cp_.get())->isEmpty());
// Distance
ensure_equals(cp_->distance(cp_.get()), 0);
diff --git a/tests/unit/geom/MultiCurveTest.cpp b/tests/unit/geom/MultiCurveTest.cpp
index 70b0ada4d..9973b5177 100644
--- a/tests/unit/geom/MultiCurveTest.cpp
+++ b/tests/unit/geom/MultiCurveTest.cpp
@@ -153,10 +153,11 @@ void object::test<3>()
// Overlay
ensure_THROW(mc_->Union(), geos::util::UnsupportedOperationException);
- ensure_THROW(mc_->Union(mc_.get()), geos::util::UnsupportedOperationException);
- ensure_THROW(mc_->difference(mc_.get()), geos::util::UnsupportedOperationException);
- ensure_THROW(mc_->intersection(mc_.get()), geos::util::UnsupportedOperationException);
- ensure_THROW(mc_->symDifference(mc_.get()), geos::util::UnsupportedOperationException);
+ // TODO: Prevent overlay from degrading CompoundCurve into MultiCurve
+ //ensure_equals_geometry(mc_->Union(mc_.get()).get(), static_cast<const Geometry*>(mc_.get()));
+ ensure(mc_->difference(mc_.get())->isEmpty());
+ //ensure_equals_geometry(mc_->intersection(mc_.get()).get(), static_cast<const Geometry*>(mc_.get()));
+ ensure(mc_->symDifference(mc_.get())->isEmpty());
// Distance
ensure_equals(mc_->distance(mc_.get()), 0);
diff --git a/tests/unit/geom/MultiSurfaceTest.cpp b/tests/unit/geom/MultiSurfaceTest.cpp
index f12fb9e38..9ee1cc244 100644
--- a/tests/unit/geom/MultiSurfaceTest.cpp
+++ b/tests/unit/geom/MultiSurfaceTest.cpp
@@ -142,10 +142,10 @@ void object::test<3>()
// Overlay
ensure_THROW(ms_->Union(), geos::util::UnsupportedOperationException);
- ensure_THROW(ms_->Union(ms_.get()), geos::util::UnsupportedOperationException);
- ensure_THROW(ms_->difference(ms_.get()), geos::util::UnsupportedOperationException);
- ensure_THROW(ms_->intersection(ms_.get()), geos::util::UnsupportedOperationException);
- ensure_THROW(ms_->symDifference(ms_.get()), geos::util::UnsupportedOperationException);
+ ensure_equals_geometry(ms_->Union(ms_.get()).get(), static_cast<const Geometry*>(ms_.get()));
+ ensure(ms_->difference(ms_.get())->isEmpty());
+ ensure_equals_geometry(ms_->intersection(ms_.get()).get(), static_cast<const Geometry*>(ms_.get()));
+ ensure(ms_->symDifference(ms_.get())->isEmpty());
// Distance
ensure_equals(ms_->distance(ms_.get()), 0);
diff --git a/tests/unit/operation/overlayng/OverlayNGTest.cpp b/tests/unit/operation/overlayng/OverlayNGTest.cpp
index 92f1d32a1..d7077e57a 100644
--- a/tests/unit/operation/overlayng/OverlayNGTest.cpp
+++ b/tests/unit/operation/overlayng/OverlayNGTest.cpp
@@ -51,20 +51,8 @@ struct test_overlayng_data {
std::unique_ptr<Geometry> geom_expected = r.read(expected);
- std::unique_ptr<Geometry> geom_result;
- if (geom_a->hasCurvedComponents() || geom_b->hasCurvedComponents()) {
- SimpleNoder noder;
- geos::algorithm::CircularArcIntersector cai;
- ArcIntersectionAdder aia(cai);
+ auto geom_result = OverlayNG::overlay(geom_a.get(), geom_b.get(), opCode, pm.get());
- noder.setArcIntersector(aia);
- OverlayNG ong(geom_a.get(), geom_b.get(), pm.get(), opCode);
- ong.setNoder(&noder);
-
- geom_result = ong.getResult();
- } else {
- geom_result = OverlayNG::overlay(geom_a.get(), geom_b.get(), opCode, pm.get());
- }
//std::string wkt_result = w.write(geom_result.get());
//std::cout << std::endl << wkt_result << std::endl;
ensure_equals_geometry_xyzm(geom_result.get(), geom_expected.get(), tol);
commit 095c72709ea5043b4a1433d2adb3611410564cd3
Author: Daniel Baston <dbaston at gmail.com>
Date: Fri Mar 27 11:31:49 2026 -0400
OverlayNG: Support curved types
diff --git a/include/geos/algorithm/PointLocator.h b/include/geos/algorithm/PointLocator.h
index be504c2b0..c249fdc42 100644
--- a/include/geos/algorithm/PointLocator.h
+++ b/include/geos/algorithm/PointLocator.h
@@ -25,6 +25,8 @@
// Forward declarations
namespace geos {
namespace geom {
+class CircularString;
+class CompoundCurve;
class CoordinateXY;
class Curve;
class Geometry;
@@ -98,6 +100,10 @@ private:
geom::Location locate(const geom::CoordinateXY& p, const geom::LineString* l);
+ geom::Location locate(const geom::CoordinateXY& p, const geom::CircularString* cs);
+
+ geom::Location locate(const geom::CoordinateXY& p, const geom::CompoundCurve* cs);
+
static geom::Location locateInPolygonRing(const geom::CoordinateXY& p, const geom::Curve* ring);
static geom::Location locate(const geom::CoordinateXY& p, const geom::Surface* poly);
diff --git a/include/geos/geom/Curve.h b/include/geos/geom/Curve.h
index 1bb16df37..8f30900f2 100644
--- a/include/geos/geom/Curve.h
+++ b/include/geos/geom/Curve.h
@@ -81,6 +81,8 @@ public:
virtual const SimpleCurve* getCurveN(std::size_t) const = 0;
+ std::unique_ptr<Curve> clone() const;
+
std::unique_ptr<Curve> reverse() const;
protected:
diff --git a/include/geos/operation/overlayng/Edge.h b/include/geos/operation/overlayng/Edge.h
index 7ec26eb16..83f67dc6a 100644
--- a/include/geos/operation/overlayng/Edge.h
+++ b/include/geos/operation/overlayng/Edge.h
@@ -184,18 +184,6 @@ private:
public:
-#if 0
- Edge()
- : aDim(OverlayLabel::DIM_UNKNOWN)
- , aDepthDelta(0)
- , aIsHole(false)
- , bDim(OverlayLabel::DIM_UNKNOWN)
- , bDepthDelta(0)
- , bIsHole(false)
- , pts(nullptr)
- {};
-#endif
-
Edge(const std::shared_ptr<const geom::CoordinateSequence>& p_pts, const EdgeSourceInfo* info, bool isCurved);
friend std::ostream& operator<<(std::ostream& os, const Edge& e);
diff --git a/include/geos/operation/overlayng/LineBuilder.h b/include/geos/operation/overlayng/LineBuilder.h
index 2abeda3be..9256f8751 100644
--- a/include/geos/operation/overlayng/LineBuilder.h
+++ b/include/geos/operation/overlayng/LineBuilder.h
@@ -77,7 +77,7 @@ private:
const geom::GeometryFactory* geometryFactory;
bool hasResultArea;
int8_t inputAreaIndex;
- std::vector<std::unique_ptr<geom::LineString>> lines;
+ std::vector<std::unique_ptr<geom::Curve>> lines;
/**
* Indicates whether intersections are allowed to produce
@@ -121,7 +121,7 @@ private:
void addResultLines();
void addResultLinesMerged();
- std::unique_ptr<geom::LineString> toLine(OverlayEdge* edge) const;
+ std::unique_ptr<geom::Curve> toLine(OverlayEdge* edge) const;
void addResultLinesForNodes();
@@ -143,7 +143,7 @@ private:
* E.g. using the direction of the majority of segments,
* or preferring the direction of the A edges.)
*/
- std::unique_ptr<geom::LineString> buildLine(OverlayEdge* node);
+ std::unique_ptr<geom::Curve> buildLine(OverlayEdge* node) const;
/*
* Finds the next edge around a node which forms
@@ -173,7 +173,7 @@ public:
LineBuilder(const LineBuilder&) = delete;
LineBuilder& operator=(const LineBuilder&) = delete;
- std::vector<std::unique_ptr<geom::LineString>> getLines();
+ std::vector<std::unique_ptr<geom::Curve>> getLines();
void setStrictMode(bool p_isStrictResultMode)
{
diff --git a/include/geos/operation/overlayng/OverlayEdge.h b/include/geos/operation/overlayng/OverlayEdge.h
index 7ad5bdf9e..85dc54fb2 100644
--- a/include/geos/operation/overlayng/OverlayEdge.h
+++ b/include/geos/operation/overlayng/OverlayEdge.h
@@ -61,6 +61,7 @@ private:
* The label must be interpreted accordingly.
*/
bool direction;
+ bool m_isCurved;
CoordinateXY dirPt;
OverlayLabel* label;
bool m_isInResultArea;
@@ -81,10 +82,12 @@ public:
OverlayEdge(const CoordinateXYZM& p_orig, const CoordinateXY& p_dirPt,
bool p_direction, OverlayLabel* p_label,
- const std::shared_ptr<const CoordinateSequence>& p_pts)
+ const std::shared_ptr<const CoordinateSequence>& p_pts,
+ bool isCurved)
: HalfEdge(p_orig)
, pts(p_pts)
, direction(p_direction)
+ , m_isCurved(isCurved)
, dirPt(p_dirPt)
, label(p_label)
, m_isInResultArea(false)
@@ -103,6 +106,11 @@ public:
return direction;
};
+ bool isCurved() const
+ {
+ return m_isCurved;
+ }
+
const CoordinateXY& directionPt() const override
{
return dirPt;
diff --git a/include/geos/operation/overlayng/OverlayEdgeRing.h b/include/geos/operation/overlayng/OverlayEdgeRing.h
index f0e94445a..c3472c15e 100644
--- a/include/geos/operation/overlayng/OverlayEdgeRing.h
+++ b/include/geos/operation/overlayng/OverlayEdgeRing.h
@@ -30,9 +30,9 @@ namespace geom {
class Coordinate;
class CoordinateXY;
class CoordinateSequence;
+class Curve;
class GeometryFactory;
-class LinearRing;
-class Polygon;
+class Surface;
}
namespace operation {
namespace overlayng {
@@ -50,8 +50,8 @@ class GEOS_DLL OverlayEdgeRing {
using CoordinateXY = geos::geom::CoordinateXY;
using CoordinateSequence = geos::geom::CoordinateSequence;
using GeometryFactory = geos::geom::GeometryFactory;
- using LinearRing = geos::geom::LinearRing;
- using Polygon = geos::geom::Polygon;
+ using Curve = geos::geom::Curve;
+ using Surface = geos::geom::Surface;
using PointOnGeometryLocator = algorithm::locate::PointOnGeometryLocator;
using IndexedPointInAreaLocator = algorithm::locate::IndexedPointInAreaLocator;
@@ -59,23 +59,22 @@ private:
// Members
OverlayEdge* startEdge;
- std::unique_ptr<LinearRing> ring;
+ std::unique_ptr<Curve> ring;
bool m_isHole;
- mutable std::unique_ptr<IndexedPointInAreaLocator> locator;
+ mutable std::unique_ptr<PointOnGeometryLocator> locator;
OverlayEdgeRing* shell;
// a list of EdgeRings which are holes in this EdgeRing
std::vector<OverlayEdgeRing*> holes;
// Methods
- void computeRingPts(OverlayEdge* start, CoordinateSequence& pts);
- void computeRing(const std::shared_ptr<CoordinateSequence> & ringPts, const GeometryFactory* geometryFactory);
+ void computeRing(OverlayEdge* start, const GeometryFactory* geometryFactory);
+ std::unique_ptr<Curve> computeRingGeometry(OverlayEdge* start, const GeometryFactory* geometryFactory) const;
/**
* Computes the list of coordinates which are contained in this ring.
* The coordinates are computed once only and cached.
* @return an array of the {@link Coordinate}s in this ring
*/
- const CoordinateSequence& getCoordinates() const;
PointOnGeometryLocator* getLocator() const;
static void closeRing(CoordinateSequence& pts);
bool contains(const OverlayEdgeRing& otherRing) const;
@@ -86,8 +85,8 @@ public:
OverlayEdgeRing(OverlayEdge* start, const GeometryFactory* geometryFactory);
- std::unique_ptr<LinearRing> getRing();
- const LinearRing* getRingPtr() const;
+ std::unique_ptr<Curve> getRing();
+ const Curve* getRingPtr() const;
const geom::Envelope& getEnvelope() const;
@@ -128,7 +127,7 @@ public:
* Computes the {@link Polygon} formed by this ring and any contained holes.
* @return the {@link Polygon} formed by this ring and its holes.
*/
- std::unique_ptr<Polygon> toPolygon(const GeometryFactory* factory);
+ std::unique_ptr<Surface> toSurface(const GeometryFactory* factory);
OverlayEdge* getEdge();
diff --git a/include/geos/operation/overlayng/OverlayGraph.h b/include/geos/operation/overlayng/OverlayGraph.h
index eda4d07ed..85a46d97a 100644
--- a/include/geos/operation/overlayng/OverlayGraph.h
+++ b/include/geos/operation/overlayng/OverlayGraph.h
@@ -70,13 +70,13 @@ private:
* Create and add HalfEdge pairs to map and vector containers,
* using local std::deque storage for objects.
*/
- OverlayEdge* createEdgePair(const std::shared_ptr<const CoordinateSequence> &pts, OverlayLabel *lbl);
+ OverlayEdge* createEdgePair(const std::shared_ptr<const CoordinateSequence> &pts, OverlayLabel *lbl, bool isCurved);
/**
* Create a single OverlayEdge in local std::deque storage, and return the
* pointer.
*/
- OverlayEdge* createOverlayEdge(const std::shared_ptr<const CoordinateSequence> &pts, OverlayLabel *lbl, bool direction);
+ OverlayEdge* createOverlayEdge(const std::shared_ptr<const CoordinateSequence> &pts, OverlayLabel *lbl, bool direction, bool isCurved);
void insert(OverlayEdge* e);
diff --git a/include/geos/operation/overlayng/OverlayMixedPoints.h b/include/geos/operation/overlayng/OverlayMixedPoints.h
index c7b681228..21907ad74 100644
--- a/include/geos/operation/overlayng/OverlayMixedPoints.h
+++ b/include/geos/operation/overlayng/OverlayMixedPoints.h
@@ -40,6 +40,9 @@ namespace locate {
class PointOnGeometryLocator;
}
}
+namespace noding {
+class Noder;
+}
}
namespace geos { // geos.
@@ -74,6 +77,7 @@ namespace overlayng { // geos.operation.overlayng
* @author Martin Davis
*/
class GEOS_DLL OverlayMixedPoints {
+ using Curve = geos::geom::Curve;
using GeometryFactory = geos::geom::GeometryFactory;
using PrecisionModel = geos::geom::PrecisionModel;
using Geometry = geos::geom::Geometry;
@@ -83,6 +87,8 @@ class GEOS_DLL OverlayMixedPoints {
using Point = geos::geom::Point;
using Polygon = geos::geom::Polygon;
using LineString = geos::geom::LineString;
+ using Noder = geos::noding::Noder;
+ using Surface = geos::geom::Surface;
using PointOnGeometryLocator = algorithm::locate::PointOnGeometryLocator;
private:
@@ -94,6 +100,7 @@ private:
const Geometry* geomNonPointInput;
const GeometryFactory* geometryFactory;
bool isPointRHS;
+ Noder* noder;
std::unique_ptr<Geometry> geomNonPoint;
int geomNonPointDim;
@@ -121,19 +128,19 @@ private:
std::unique_ptr<Geometry> copyNonPoint() const;
- std::unique_ptr<CoordinateSequence> extractCoordinates(const Geometry* points, const PrecisionModel* pm) const;
+ static std::unique_ptr<CoordinateSequence> extractCoordinates(const Geometry* points, const PrecisionModel* pm);
- std::vector<std::unique_ptr<Polygon>> extractPolygons(const Geometry* geom) const;
+ static std::vector<std::unique_ptr<Surface>> extractPolygons(const Geometry* geom);
- std::vector<std::unique_ptr<LineString>> extractLines(const Geometry* geom) const;
+ static std::vector<std::unique_ptr<Curve>> extractLines(const Geometry* geom);
public:
- OverlayMixedPoints(int p_opCode, const Geometry* geom0, const Geometry* geom1, const PrecisionModel* p_pm);
+ OverlayMixedPoints(int p_opCode, const Geometry* geom0, const Geometry* geom1, const PrecisionModel* p_pm, Noder* p_noder);
- static std::unique_ptr<Geometry> overlay(int opCode, const Geometry* geom0, const Geometry* geom1, const PrecisionModel* pm);
+ static std::unique_ptr<Geometry> overlay(int opCode, const Geometry* geom0, const Geometry* geom1, const PrecisionModel* pm, Noder* p_noder);
std::unique_ptr<Geometry> getResult();
diff --git a/include/geos/operation/overlayng/OverlayUtil.h b/include/geos/operation/overlayng/OverlayUtil.h
index 6bff9c358..25fa99312 100644
--- a/include/geos/operation/overlayng/OverlayUtil.h
+++ b/include/geos/operation/overlayng/OverlayUtil.h
@@ -14,10 +14,10 @@
#pragma once
-#include <geos/geom/Point.h>
-#include <geos/geom/Polygon.h>
-#include <geos/geom/LineString.h>
+#include <geos/geom/Curve.h>
#include <geos/geom/Geometry.h>
+#include <geos/geom/Point.h>
+#include <geos/geom/Surface.h>
#include <geos/export.h>
@@ -55,10 +55,10 @@ class GEOS_DLL OverlayUtil {
using Geometry = geos::geom::Geometry;
using Coordinate = geos::geom::Coordinate;
using CoordinateSequence = geos::geom::CoordinateSequence;
+ using Curve = geos::geom::Curve;
using Envelope = geos::geom::Envelope;
using Point = geos::geom::Point;
- using LineString = geos::geom::LineString;
- using Polygon = geos::geom::Polygon;
+ using Surface = geos::geom::Surface;
using GeometryFactory = geos::geom::GeometryFactory;
using PrecisionModel = geos::geom::PrecisionModel;
@@ -172,8 +172,8 @@ public:
* Creates an overlay result geometry for homogeneous or mixed components.
*/
static std::unique_ptr<Geometry> createResultGeometry(
- std::vector<std::unique_ptr<Polygon>>& resultPolyList,
- std::vector<std::unique_ptr<LineString>>& resultLineList,
+ std::vector<std::unique_ptr<Surface>>& resultPolyList,
+ std::vector<std::unique_ptr<Curve>>& resultLineList,
std::vector<std::unique_ptr<Point>>& resultPointList,
const GeometryFactory* geometryFactory);
diff --git a/include/geos/operation/overlayng/PolygonBuilder.h b/include/geos/operation/overlayng/PolygonBuilder.h
index 9e820fc14..f30b43226 100644
--- a/include/geos/operation/overlayng/PolygonBuilder.h
+++ b/include/geos/operation/overlayng/PolygonBuilder.h
@@ -29,7 +29,7 @@ namespace geos {
namespace geom {
class GeometryFactory;
class Geometry;
-class Polygon;
+class Surface;
}
namespace operation {
namespace overlayng {
@@ -55,7 +55,7 @@ private:
// Storage
std::vector<std::unique_ptr<OverlayEdgeRing>> vecOER;
- std::vector<std::unique_ptr<geom::Polygon>> computePolygons(const std::vector<OverlayEdgeRing*>& shellList) const;
+ std::vector<std::unique_ptr<geom::Surface>> computePolygons(const std::vector<OverlayEdgeRing*>& shellList) const;
void buildRings(const std::vector<OverlayEdge*>& resultAreaEdges);
@@ -138,7 +138,7 @@ public:
PolygonBuilder& operator=(const PolygonBuilder&) = delete;
// Methods
- std::vector<std::unique_ptr<geom::Polygon>> getPolygons() const;
+ std::vector<std::unique_ptr<geom::Surface>> getPolygons() const;
std::vector<OverlayEdgeRing*> getShellRings() const;
diff --git a/src/algorithm/PointLocator.cpp b/src/algorithm/PointLocator.cpp
index 5545d7655..abcee3e12 100644
--- a/src/algorithm/PointLocator.cpp
+++ b/src/algorithm/PointLocator.cpp
@@ -19,6 +19,8 @@
#include <geos/algorithm/PointLocator.h>
#include <geos/algorithm/PointLocation.h>
+#include <geos/geom/CircularString.h>
+#include <geos/geom/CompoundCurve.h>
#include <geos/geom/Geometry.h>
#include <geos/geom/Point.h>
#include <geos/geom/LineString.h>
@@ -31,7 +33,6 @@
#include <geos/util/UnsupportedOperationException.h>
#include <cassert>
-#include <typeinfo>
using namespace geos::geom;
@@ -97,6 +98,16 @@ PointLocator::computeLocation(const CoordinateXY& p, const Geometry* geom)
updateLocationInfo(locate(p, ls));
break;
}
+ case GEOS_CIRCULARSTRING: {
+ const CircularString* cs = static_cast<const CircularString*>(geom);
+ updateLocationInfo(locate(p, cs));
+ break;
+ }
+ case GEOS_COMPOUNDCURVE: {
+ const CompoundCurve* cc = static_cast<const CompoundCurve*>(geom);
+ updateLocationInfo(locate(p, cc));
+ break;
+ }
case GEOS_POLYGON: {
const Polygon* po = static_cast<const Polygon*>(geom);
updateLocationInfo(locate(p, po));
@@ -183,6 +194,73 @@ PointLocator::locate(const CoordinateXY& p, const LineString* l)
return Location::EXTERIOR;
}
+/* private */
+Location
+PointLocator::locate(const CoordinateXY& p, const CircularString* cs)
+{
+ if(!cs->getEnvelopeInternal()->intersects(p)) {
+ return Location::EXTERIOR;
+ }
+
+ if (cs->isClosed()) {
+ const CoordinateSequence* seq = cs->getCoordinatesRO();
+ if (p.equals2D(seq->front<CoordinateXY>())) {
+ return Location::BOUNDARY;
+ }
+ if (p.equals2D(seq->back<CoordinateXY>())) {
+ return Location::BOUNDARY;
+ }
+ }
+
+ for (const auto& arc : cs->getArcs()) {
+ if (arc.containsPoint(p)) {
+ return Location::INTERIOR;
+ }
+ }
+
+ return Location::EXTERIOR;
+}
+
+/* private */
+Location
+PointLocator::locate(const CoordinateXY& p, const CompoundCurve* cc)
+{
+ if(!cc->getEnvelopeInternal()->intersects(p)) {
+ return Location::EXTERIOR;
+ }
+
+ if (cc->isClosed()) {
+ const SimpleCurve* firstCurve = cc->getCurveN(0);
+ if (p.equals2D(firstCurve->getCoordinatesRO()->front<CoordinateXY>())) {
+ return Location::BOUNDARY;
+ }
+
+ const SimpleCurve* lastCurve = cc->getCurveN(cc->getNumCurves() - 1);
+ if (p.equals2D(lastCurve->getCoordinatesRO()->back<CoordinateXY>())) {
+ return Location::BOUNDARY;
+ }
+ }
+
+ for (std::size_t i = 0; i < cc->getNumCurves(); i++) {
+ const SimpleCurve* curve = cc->getCurveN(i);
+ if (curve->getGeometryTypeId() == GEOS_CIRCULARSTRING) {
+ const CircularString* cs = static_cast<const CircularString*>(curve);
+ for (const auto& arc : cs->getArcs()) {
+ if (arc.containsPoint(p)) {
+ return Location::INTERIOR;
+ }
+ }
+ } else {
+ const CoordinateSequence* seq = curve->getCoordinatesRO();
+ if(PointLocation::isOnLine(p, seq)) {
+ return Location::INTERIOR;
+ }
+ }
+ }
+
+ return Location::EXTERIOR;
+}
+
/* private */
Location
PointLocator::locateInPolygonRing(const CoordinateXY& p, const Curve* ring)
diff --git a/src/geom/Curve.cpp b/src/geom/Curve.cpp
index a55678267..26b7dc533 100644
--- a/src/geom/Curve.cpp
+++ b/src/geom/Curve.cpp
@@ -48,6 +48,12 @@ Curve::apply_rw(GeometryFilter* filter)
filter->filter_rw(this);
}
+std::unique_ptr<Curve>
+Curve::clone() const
+{
+ return std::unique_ptr<Curve>(static_cast<Curve*>(Geometry::clone().release()));
+}
+
bool
Curve::isRing() const
{
diff --git a/src/geom/GeometryFactory.cpp b/src/geom/GeometryFactory.cpp
index f53012a4c..cc5e97870 100644
--- a/src/geom/GeometryFactory.cpp
+++ b/src/geom/GeometryFactory.cpp
@@ -642,6 +642,13 @@ const
return std::unique_ptr<CurvePolygon>(new CurvePolygon(std::move(shell), std::move(holes), *this));
}
+static std::unique_ptr<LineString>
+linearRingToLineString(std::unique_ptr<Curve> linearRing)
+{
+ const std::shared_ptr<const CoordinateSequence> pts = detail::down_cast<const SimpleCurve*>(linearRing.get())->getSharedCoordinates();
+ return linearRing->getFactory()->createLineString(pts);
+}
+
/* public */
std::unique_ptr<Surface>
GeometryFactory::createSurface(std::unique_ptr<Curve> && shell)
@@ -652,6 +659,10 @@ const
return createPolygon(std::move(shellLR));
}
+ if (shell->getGeometryTypeId() == GEOS_LINEARRING) {
+ shell = linearRingToLineString(std::move(shell));
+ }
+
return createCurvePolygon(std::move(shell));
}
@@ -675,6 +686,17 @@ const
return createPolygon(std::move(shellLR), std::move(holesLR));
}
+ // Purge LinearRings
+ if (shell->getGeometryTypeId() == GEOS_LINEARRING) {
+ shell = linearRingToLineString(std::move(shell));
+ }
+
+ for (std::size_t i = 0; i < holes.size(); i++) {
+ if (holes[i]->getGeometryTypeId() == GEOS_LINEARRING) {
+ holes[i] = linearRingToLineString(std::move(holes[i]));
+ }
+ }
+
return createCurvePolygon(std::move(shell), std::move(holes));
}
diff --git a/src/noding/snap/SnappingNoder.cpp b/src/noding/snap/SnappingNoder.cpp
index f4848d1fc..e2f0643d6 100644
--- a/src/noding/snap/SnappingNoder.cpp
+++ b/src/noding/snap/SnappingNoder.cpp
@@ -111,9 +111,8 @@ SnappingNoder::snap(const CoordinateSequence *cs)
snapCoords->reserve(cs->size());
cs->forEach<Coordinate>([&snapCoords, this](const Coordinate& origPt) {
- const Coordinate& pt = snapIndex.snap(origPt);
+ const Coordinate& pt = origPt.isValid() ? snapIndex.snap(origPt) : origPt;
snapCoords->add(pt, false); // Remove repeated points
-
});
return snapCoords;
}
diff --git a/src/operation/overlayng/IndexedPointOnLineLocator.cpp b/src/operation/overlayng/IndexedPointOnLineLocator.cpp
index 95a7e12da..b521145cb 100644
--- a/src/operation/overlayng/IndexedPointOnLineLocator.cpp
+++ b/src/operation/overlayng/IndexedPointOnLineLocator.cpp
@@ -29,7 +29,7 @@ namespace overlayng { // geos.operation.overlayng
/*public*/
geom::Location
IndexedPointOnLineLocator::locate(const geom::CoordinateXY* p) {
- // TODO: optimize this with a segment index
+ // TODO: optimize this with a segment/arc index
algorithm::PointLocator locator;
return locator.locate(*p, &inputGeom);
}
diff --git a/src/operation/overlayng/InputGeometry.cpp b/src/operation/overlayng/InputGeometry.cpp
index 6c3c5a2ba..320d887f0 100644
--- a/src/operation/overlayng/InputGeometry.cpp
+++ b/src/operation/overlayng/InputGeometry.cpp
@@ -13,6 +13,7 @@
**********************************************************************/
#include <geos/operation/overlayng/InputGeometry.h>
+#include <geos/algorithm/locate/SimplePointInAreaLocator.h>
namespace geos { // geos
namespace operation { // geos.operation
@@ -21,6 +22,7 @@ namespace overlayng { // geos.operation.overlayng
using geos::geom::Location;
using geos::geom::Geometry;
using geos::geom::Envelope;
+using geos::algorithm::locate::SimplePointInAreaLocator;
using geos::algorithm::locate::IndexedPointInAreaLocator;
using geos::algorithm::locate::PointOnGeometryLocator;
@@ -166,13 +168,23 @@ PointOnGeometryLocator*
InputGeometry::getLocator(uint8_t geomIndex)
{
if (geomIndex == 0) {
- if (ptLocatorA == nullptr)
- ptLocatorA.reset(new IndexedPointInAreaLocator(*getGeometry(geomIndex)));
+ if (ptLocatorA == nullptr) {
+ if (geom[0]->hasCurvedComponents()) {
+ ptLocatorA = std::make_unique<SimplePointInAreaLocator>(*getGeometry(geomIndex));
+ } else {
+ ptLocatorA = std::make_unique<IndexedPointInAreaLocator>(*getGeometry(geomIndex));
+ }
+ }
return ptLocatorA.get();
}
else {
- if (ptLocatorB == nullptr)
- ptLocatorB.reset(new IndexedPointInAreaLocator(*getGeometry(geomIndex)));
+ if (ptLocatorB == nullptr) {
+ if (geom[1]->hasCurvedComponents()) {
+ ptLocatorB = std::make_unique<SimplePointInAreaLocator>(*getGeometry(geomIndex));
+ } else {
+ ptLocatorB = std::make_unique<IndexedPointInAreaLocator>(*getGeometry(geomIndex));
+ }
+ }
return ptLocatorB.get();
}
}
diff --git a/src/operation/overlayng/LineBuilder.cpp b/src/operation/overlayng/LineBuilder.cpp
index a676b36e5..4d1a25877 100644
--- a/src/operation/overlayng/LineBuilder.cpp
+++ b/src/operation/overlayng/LineBuilder.cpp
@@ -17,8 +17,10 @@
#include <geos/operation/overlayng/OverlayEdge.h>
#include <geos/operation/overlayng/OverlayLabel.h>
#include <geos/operation/overlayng/OverlayNG.h>
+#include <geos/geom/CircularString.h>
#include <geos/geom/GeometryFactory.h>
#include <geos/geom/CoordinateSequence.h>
+#include <geos/geom/util/CurveBuilder.h>
@@ -29,7 +31,7 @@ namespace overlayng { // geos.operation.overlayng
using namespace geos::geom;
/*public*/
-std::vector<std::unique_ptr<LineString>>
+std::vector<std::unique_ptr<Curve>>
LineBuilder::getLines()
{
markResultLines();
@@ -174,19 +176,22 @@ LineBuilder::addResultLines()
}
}
-std::unique_ptr<LineString>
+std::unique_ptr<Curve>
LineBuilder::toLine(OverlayEdge* edge) const
{
// bool isForward = edge->isForward();
const auto* edgePts = edge->getCoordinatesRO();
- std::unique_ptr<CoordinateSequence> pts(new CoordinateSequence(0u, edgePts->hasZ(), edgePts->hasM()));
+ auto pts = std::make_unique<CoordinateSequence>(0u, edgePts->hasZ(), edgePts->hasM());
pts->reserve(edgePts->size());
pts->add(edge->orig(), false);
edge->addCoordinates(pts.get());
assert(pts->size() == edgePts->size());
+ if (edge->isCurved()) {
+ return geometryFactory->createCircularString(std::move(pts));
+ }
return geometryFactory->createLineString(std::move(pts));
}
@@ -217,7 +222,7 @@ LineBuilder::addResultLinesForNodes()
* This will find all lines originating at nodes
*/
if (degreeOfLines(edge) != 2) {
- std::unique_ptr<LineString> line = buildLine(edge);
+ auto line = buildLine(edge);
lines.push_back(std::move(line));
}
}
@@ -241,20 +246,24 @@ LineBuilder::addResultLinesRings()
}
/*private*/
-std::unique_ptr<LineString>
-LineBuilder::buildLine(OverlayEdge* node)
+std::unique_ptr<Curve>
+LineBuilder::buildLine(OverlayEdge* node) const
{
+ const bool constructZ = node->getCoordinatesRO()->hasZ();
+ const bool constructM = node->getCoordinatesRO()->hasM();
+ geom::util::CurveBuilder cb(*geometryFactory, constructZ, constructM);
+
// assert: edgeStart degree = 1
// assert: edgeStart direction = forward
- std::unique_ptr<CoordinateSequence> pts(new CoordinateSequence());
- pts->add(node->orig(), false);
+ cb.add(*node->getCoordinatesRO(), node->isCurved());
bool isNodeForward = node->isForward();
OverlayEdge* e = node;
do {
e->markVisitedBoth();
- e->addCoordinates(pts.get());
+ CoordinateSequence& pts = cb.getSeq(e->isCurved());
+ e->addCoordinates(&pts);
// end line if next vertex is a node
if (degreeOfLines(e->symOE()) != 2) {
@@ -264,12 +273,14 @@ LineBuilder::buildLine(OverlayEdge* node)
// e will be nullptr if next edge has been visited, which indicates a ring
}
while (e != nullptr);
+
// reverse coordinates before constructing
+ auto geom = cb.getGeometry();
if(!isNodeForward) {
- pts->reverse();
+ return geom->reverse();
}
- return geometryFactory->createLineString(std::move(pts));
+ return geom;
}
/*private*/
diff --git a/src/operation/overlayng/OverlayEdgeRing.cpp b/src/operation/overlayng/OverlayEdgeRing.cpp
index 1587f9e0c..12f4c8e92 100644
--- a/src/operation/overlayng/OverlayEdgeRing.cpp
+++ b/src/operation/overlayng/OverlayEdgeRing.cpp
@@ -17,11 +17,18 @@
#include <geos/geom/Location.h>
#include <geos/geom/Envelope.h>
#include <geos/geom/Coordinate.h>
+#include <geos/geom/CoordinateFilter.h>
#include <geos/geom/CoordinateSequence.h>
+#include <geos/geom/CircularString.h>
+#include <geos/geom/CompoundCurve.h>
+#include <geos/geom/CurvePolygon.h>
#include <geos/geom/GeometryFactory.h>
+#include <geos/geom/util/CurveBuilder.h>
#include <geos/util/TopologyException.h>
+#include <geos/algorithm/PointLocation.h>
#include <geos/algorithm/locate/PointOnGeometryLocator.h>
#include <geos/algorithm/locate/IndexedPointInAreaLocator.h>
+#include <geos/algorithm/locate/SimplePointInAreaLocator.h>
#include <geos/algorithm/Orientation.h>
#include <geos/operation/polygonize/EdgeRing.h>
@@ -32,8 +39,7 @@ namespace overlayng { // geos.operation.overlayng
using namespace geos::geom;
using geos::operation::polygonize::EdgeRing;
using geos::algorithm::locate::PointOnGeometryLocator;
-
-
+using geos::algorithm::locate::SimplePointInAreaLocator;
OverlayEdgeRing::OverlayEdgeRing(OverlayEdge* start, const GeometryFactory* geometryFactory)
: startEdge(start)
@@ -42,19 +48,17 @@ OverlayEdgeRing::OverlayEdgeRing(OverlayEdge* start, const GeometryFactory* geom
, locator(nullptr)
, shell(nullptr)
{
- auto ringPts = std::make_shared<CoordinateSequence>(0u, start->getCoordinatesRO()->hasZ(), start->getCoordinatesRO()->hasM());
- computeRingPts(start, *ringPts);
- computeRing(ringPts, geometryFactory);
+ computeRing(start, geometryFactory);
}
/*public*/
-std::unique_ptr<LinearRing>
+std::unique_ptr<Curve>
OverlayEdgeRing::getRing()
{
return std::move(ring);
}
-const LinearRing*
+const Curve*
OverlayEdgeRing::getRingPtr() const
{
return ring.get();
@@ -133,46 +137,49 @@ OverlayEdgeRing::closeRing(CoordinateSequence& pts)
}
/*private*/
-void
-OverlayEdgeRing::computeRingPts(OverlayEdge* start, CoordinateSequence& pts)
+std::unique_ptr<Curve>
+OverlayEdgeRing::computeRingGeometry(OverlayEdge* start, const GeometryFactory* gfact) const
{
OverlayEdge* edge = start;
+
+ const bool hasZ = start->getCoordinatesRO()->hasZ();
+ const bool hasM = start->getCoordinatesRO()->hasM();
+
+ geom::util::CurveBuilder cb(*gfact, hasZ, hasM);
+
do {
if (edge->getEdgeRing() == this)
- throw util::TopologyException("Edge visited twice during ring-building", edge->getCoordinate());
+ throw geos::util::TopologyException("Edge visited twice during ring-building", edge->getCoordinate());
// only valid for polygonal output
+ CoordinateSequence& pts = cb.getSeq(edge->isCurved());
edge->addCoordinates(&pts);
edge->setEdgeRing(this);
+
if (edge->nextResult() == nullptr)
- throw util::TopologyException("Found null edge in ring", edge->dest());
+ throw geos::util::TopologyException("Found null edge in ring", edge->dest());
edge = edge->nextResult();
}
while (edge != start);
- closeRing(pts);
+
+ cb.closeRing();
+ return cb.getGeometry();
}
/*private*/
void
-OverlayEdgeRing::computeRing(const std::shared_ptr<CoordinateSequence> & p_ringPts, const GeometryFactory* geometryFactory)
+OverlayEdgeRing::computeRing(OverlayEdge* start, const GeometryFactory* geometryFactory)
{
if (ring != nullptr) return; // don't compute more than once
- ring = geometryFactory->createLinearRing(p_ringPts);
- m_isHole = algorithm::Orientation::isCCW(ring->getCoordinatesRO());
-}
+ ring = computeRingGeometry(start, geometryFactory);
-/**
-* Computes the list of coordinates which are contained in this ring.
-* The coordinates are computed once only and cached.
-*
-* @return an array of the {@link Coordinate}s in this ring
-*/
-/*private*/
-const CoordinateSequence&
-OverlayEdgeRing::getCoordinates() const
-{
- return *ring->getCoordinatesRO();
+ if (ring->getGeometryTypeId() == GEOS_COMPOUNDCURVE) {
+ const auto seq = ring->getCoordinates();
+ m_isHole = algorithm::Orientation::isCCW(seq.get());
+ } else {
+ m_isHole = algorithm::Orientation::isCCW(static_cast<const SimpleCurve*>(ring.get())->getCoordinatesRO());
+ }
}
/**
@@ -209,12 +216,33 @@ OverlayEdgeRing::findEdgeRingContaining(const std::vector<OverlayEdgeRing*>& erL
return minContainingRing;
}
+// Adapter class to check whether a point lies within a ring.
+// Unlike IndexedPointInAreaLocator, SimplePointInAreaLocator does not treat
+// closed rings as areas.
+class PointInCurvedRingLocator : public algorithm::locate::PointOnGeometryLocator {
+public:
+ explicit PointInCurvedRingLocator(const Curve& ring) : m_ring(ring) {}
+
+ geom::Location locate(const geom::CoordinateXY* p) override {
+ return algorithm::PointLocation::locateInRing(*p, m_ring);
+ }
+
+private:
+ const Curve& m_ring;
+};
+
/*private*/
PointOnGeometryLocator*
OverlayEdgeRing::getLocator() const
{
if (locator == nullptr) {
- locator.reset(new IndexedPointInAreaLocator(*(getRingPtr())));
+ const Curve* p_ring = getRingPtr();
+ if (p_ring->hasCurvedComponents()) {
+ locator = detail::make_unique<PointInCurvedRingLocator>(*p_ring);
+ } else {
+ locator = detail::make_unique<IndexedPointInAreaLocator>(*p_ring);
+ }
+
}
return locator.get();
}
@@ -252,24 +280,51 @@ OverlayEdgeRing::contains(const OverlayEdgeRing& otherRing) const {
bool
OverlayEdgeRing::isPointInOrOut(const OverlayEdgeRing& otherRing) const {
// in most cases only one or two points will be checked
- for (const CoordinateXY& pt : otherRing.getCoordinates().items<CoordinateXY>()) {
- geom::Location loc = locate(pt);
- if (loc == geom::Location::INTERIOR) {
- return true;
+
+ struct PointTester : public geom::CoordinateFilter {
+
+ public:
+ explicit PointTester(const OverlayEdgeRing* p_ring) : m_er(p_ring), m_result(false) {}
+
+ void filter_ro(const CoordinateXY* pt) override {
+ Location loc = m_er->locate(*pt);
+
+ if (loc == geom::Location::INTERIOR) {
+ m_done = true;
+ m_result = true;
+ } else if (loc == geom::Location::EXTERIOR) {
+ m_done = true;
+ m_result = false;
+ }
+
+ // pt is on BOUNDARY, so keep checking for a determining location
}
- if (loc == geom::Location::EXTERIOR) {
- return false;
+
+ bool isDone() const override {
+ return m_done;
}
- // pt is on BOUNDARY, so keep checking for a determining location
- }
- return false;
+
+ bool getResult() const {
+ return m_result;
+ }
+
+ private:
+ const OverlayEdgeRing* m_er;
+ bool m_done{false};
+ bool m_result;
+ };
+
+ PointTester tester(this);
+ otherRing.getRingPtr()->apply_ro(&tester);
+
+ return tester.getResult();
}
/*public*/
const Coordinate&
OverlayEdgeRing::getCoordinate() const
{
- return ring->getCoordinatesRO()->getAt(0);
+ return *detail::down_cast<const Coordinate*>(ring->getCoordinate());
}
/**
@@ -277,20 +332,19 @@ OverlayEdgeRing::getCoordinate() const
* @return the {@link Polygon} formed by this ring and its holes.
*/
/*public*/
-std::unique_ptr<Polygon>
-OverlayEdgeRing::toPolygon(const GeometryFactory* factory)
+std::unique_ptr<Surface>
+OverlayEdgeRing::toSurface(const GeometryFactory* factory)
{
if (holes.empty()) {
- return factory->createPolygon(std::move(ring));
- } else {
- std::vector<std::unique_ptr<LinearRing>> holeLR(holes.size());
-
- for (std::size_t i = 0; i < holes.size(); i++) {
- holeLR[i] = holes[i]->getRing();
- }
-
- return factory->createPolygon(std::move(ring), std::move(holeLR));
+ return factory->createSurface(std::move(ring));
}
+
+ std::vector<std::unique_ptr<Curve>> holeCurves(holes.size());
+ for (std::size_t i = 0; i < holes.size(); i++) {
+ holeCurves[i] = holes[i]->getRing();
+ }
+
+ return factory->createSurface(std::move(ring), std::move(holeCurves));
}
/*public*/
diff --git a/src/operation/overlayng/OverlayGraph.cpp b/src/operation/overlayng/OverlayGraph.cpp
index a58983e21..7fd39fefd 100644
--- a/src/operation/overlayng/OverlayGraph.cpp
+++ b/src/operation/overlayng/OverlayGraph.cpp
@@ -14,9 +14,11 @@
#include <geos/operation/overlayng/OverlayGraph.h>
+#include <geos/algorithm/CircularArcs.h>
#include <geos/operation/overlayng/Edge.h>
#include <geos/geom/Coordinate.h>
#include <geos/geom/CoordinateSequence.h>
+#include <geos/geom/CircularArc.h>
#ifndef GEOS_DEBUG
#define GEOS_DEBUG 0
@@ -84,7 +86,7 @@ OverlayGraph::addEdge(Edge* edge)
{
// CoordinateSequence* pts = = edge->getCoordinates().release();
std::shared_ptr<const CoordinateSequence> pts = edge->releaseCoordinates();
- OverlayEdge* e = createEdgePair(pts, createOverlayLabel(edge));
+ OverlayEdge* e = createEdgePair(pts, createOverlayLabel(edge), edge->isCurved());
#if GEOS_DEBUG
std::cerr << "added edge: " << *e << std::endl;
#endif
@@ -95,32 +97,52 @@ OverlayGraph::addEdge(Edge* edge)
/*private*/
OverlayEdge*
-OverlayGraph::createEdgePair(const std::shared_ptr<const CoordinateSequence>& pts, OverlayLabel *lbl)
+OverlayGraph::createEdgePair(const std::shared_ptr<const CoordinateSequence>& pts, OverlayLabel *lbl, bool isCurved)
{
- OverlayEdge* e0 = createOverlayEdge(pts, lbl, true);
- OverlayEdge* e1 = createOverlayEdge(pts, lbl, false);
+ OverlayEdge* e0 = createOverlayEdge(pts, lbl, true, isCurved);
+ OverlayEdge* e1 = createOverlayEdge(pts, lbl, false, isCurved);
e0->link(e1);
return e0;
}
+static CoordinateXY
+getDirectionPoint(const CoordinateSequence& pts, bool forward, bool isCurved)
+{
+ if (isCurved) {
+ assert(pts.size() >= 3);
+ if (forward) {
+ CircularArc arc(pts, 0);
+ return arc.getDirectionPoint();
+ } else {
+ CircularArc arc(pts, pts.size() - 3);
+ return algorithm::CircularArcs::getDirectionPoint(arc.getCenter(), arc.getRadius(), arc.theta2(), !arc.isCCW());
+ }
+ }
+
+ assert(pts.size() >= 2);
+ if (forward) {
+ return pts.getAt<CoordinateXY>(1);
+ }
+
+ return pts.getAt<CoordinateXY>(pts.size() - 2);
+}
+
/*private*/
OverlayEdge*
-OverlayGraph::createOverlayEdge(const std::shared_ptr<const CoordinateSequence>& pts, OverlayLabel* lbl, bool direction)
+OverlayGraph::createOverlayEdge(const std::shared_ptr<const CoordinateSequence>& pts, OverlayLabel* lbl, bool direction, bool isCurved)
{
+ assert(pts->size() >= 2);
+
CoordinateXYZM origin;
- CoordinateXYZM dirPt;
+ const CoordinateXY dirPt = getDirectionPoint(*pts, direction, isCurved);
if (direction) {
pts->getAt(0, origin);
- pts->getAt(1, dirPt);
}
else {
- assert(pts->size() > 0);
- std::size_t ilast = pts->size() - 1;
- pts->getAt(ilast, origin);
- pts->getAt(ilast-1, dirPt);
+ pts->getAt(pts->size() - 1, origin);
}
- ovEdgeQue.emplace_back(origin, dirPt, direction, lbl, pts);
+ ovEdgeQue.emplace_back(origin, dirPt, direction, lbl, pts, isCurved);
OverlayEdge& ove = ovEdgeQue.back();
return &ove;
}
diff --git a/src/operation/overlayng/OverlayMixedPoints.cpp b/src/operation/overlayng/OverlayMixedPoints.cpp
index 62492d6b7..d5a0b86e2 100644
--- a/src/operation/overlayng/OverlayMixedPoints.cpp
+++ b/src/operation/overlayng/OverlayMixedPoints.cpp
@@ -16,6 +16,7 @@
#include <geos/algorithm/locate/IndexedPointInAreaLocator.h>
#include <geos/algorithm/locate/PointOnGeometryLocator.h>
+#include <geos/algorithm/locate/SimplePointInAreaLocator.h>
#include <geos/geom/CoordinateSequence.h>
#include <geos/geom/Coordinate.h>
#include <geos/geom/CoordinateFilter.h>
@@ -38,6 +39,7 @@ namespace overlayng { // geos.operation.overlayng
using namespace geos::geom;
using geos::algorithm::locate::PointOnGeometryLocator;
using geos::algorithm::locate::IndexedPointInAreaLocator;
+using geos::algorithm::locate::SimplePointInAreaLocator;
/**
* @brief Extracts and rounds coordinates from a geometry
@@ -76,10 +78,11 @@ private:
};
/*public*/
-OverlayMixedPoints::OverlayMixedPoints(int p_opCode, const Geometry* geom0, const Geometry* geom1, const PrecisionModel* p_pm)
+OverlayMixedPoints::OverlayMixedPoints(int p_opCode, const Geometry* geom0, const Geometry* geom1, const PrecisionModel* p_pm, Noder* p_noder)
: opCode(p_opCode)
, pm(p_pm ? p_pm : geom0->getPrecisionModel())
, geometryFactory(geom0->getFactory())
+ , noder(p_noder)
, resultDim(OverlayUtil::resultDimension(opCode, geom0->getDimension(), geom1->getDimension()))
{
// name the dimensional geometries
@@ -97,9 +100,9 @@ OverlayMixedPoints::OverlayMixedPoints(int p_opCode, const Geometry* geom0, cons
/*public static*/
std::unique_ptr<Geometry>
-OverlayMixedPoints::overlay(int opCode, const Geometry* geom0, const Geometry* geom1, const PrecisionModel* pm)
+OverlayMixedPoints::overlay(int opCode, const Geometry* geom0, const Geometry* geom1, const PrecisionModel* pm, Noder* p_noder)
{
- OverlayMixedPoints overlay(opCode, geom0, geom1, pm);
+ OverlayMixedPoints overlay(opCode, geom0, geom1, pm, p_noder);
return overlay.getResult();
}
@@ -137,8 +140,11 @@ std::unique_ptr<PointOnGeometryLocator>
OverlayMixedPoints::createLocator(const Geometry* p_geomNonPoint)
{
if (geomNonPointDim == 2) {
- std::unique_ptr<PointOnGeometryLocator> ipial(new IndexedPointInAreaLocator(*p_geomNonPoint));
- return ipial;
+ if (p_geomNonPoint->hasCurvedComponents()) {
+ return std::make_unique<SimplePointInAreaLocator>(*p_geomNonPoint);
+ } else {
+ return std::make_unique<IndexedPointInAreaLocator>(*p_geomNonPoint);
+ }
}
else {
std::unique_ptr<PointOnGeometryLocator> ipoll(new IndexedPointOnLineLocator(*p_geomNonPoint));
@@ -159,7 +165,7 @@ OverlayMixedPoints::prepareNonPoint(const Geometry* geomInput)
return geomInput->clone();
}
// Node and round the non-point geometry for output
- return OverlayNG::geomunion(geomInput, pm);
+ return OverlayNG::geomunion(geomInput, pm, noder);
}
/*private*/
@@ -175,11 +181,11 @@ std::unique_ptr<Geometry>
OverlayMixedPoints::computeUnion(const CoordinateSequence* coords)
{
std::vector<std::unique_ptr<Point>> resultPointList = findPoints(false, coords);
- std::vector<std::unique_ptr<LineString>> resultLineList;
+ std::vector<std::unique_ptr<Curve>> resultLineList;
if (geomNonPointDim == 1) {
resultLineList = extractLines(geomNonPoint.get());
}
- std::vector<std::unique_ptr<Polygon>> resultPolyList;
+ std::vector<std::unique_ptr<Surface>> resultPolyList;
if (geomNonPointDim == 2) {
resultPolyList = extractPolygons(geomNonPoint.get());
}
@@ -274,7 +280,7 @@ OverlayMixedPoints::hasLocation(bool isCovered, const CoordinateXY& coord) const
/*private*/
std::unique_ptr<CoordinateSequence>
-OverlayMixedPoints::extractCoordinates(const Geometry* points, const PrecisionModel* p_pm) const
+OverlayMixedPoints::extractCoordinates(const Geometry* points, const PrecisionModel* p_pm)
{
auto coords = detail::make_unique<CoordinateSequence>(0u, points->hasZ(), points->hasM());
coords->reserve(points->getNumPoints());
@@ -285,12 +291,12 @@ OverlayMixedPoints::extractCoordinates(const Geometry* points, const PrecisionMo
}
/*private*/
-std::vector<std::unique_ptr<Polygon>>
-OverlayMixedPoints::extractPolygons(const Geometry* geom) const
+std::vector<std::unique_ptr<Surface>>
+OverlayMixedPoints::extractPolygons(const Geometry* geom)
{
- std::vector<std::unique_ptr<Polygon>> list;
+ std::vector<std::unique_ptr<Surface>> list;
for (std::size_t i = 0; i < geom->getNumGeometries(); i++) {
- const Polygon* poly = static_cast<const Polygon*>(geom->getGeometryN(i));
+ const Surface* poly = static_cast<const Surface*>(geom->getGeometryN(i));
if(!poly->isEmpty()) {
list.emplace_back(poly->clone());
}
@@ -299,12 +305,12 @@ OverlayMixedPoints::extractPolygons(const Geometry* geom) const
}
/*private*/
-std::vector<std::unique_ptr<LineString>>
-OverlayMixedPoints::extractLines(const Geometry* geom) const
+std::vector<std::unique_ptr<Curve>>
+OverlayMixedPoints::extractLines(const Geometry* geom)
{
- std::vector<std::unique_ptr<LineString>> list;
+ std::vector<std::unique_ptr<Curve>> list;
for (std::size_t i = 0; i < geom->getNumGeometries(); i++) {
- const LineString* line = static_cast<const LineString*>(geom->getGeometryN(i));
+ const Curve* line = static_cast<const Curve*>(geom->getGeometryN(i));
if (! line->isEmpty()) {
list.emplace_back(line->clone());
}
diff --git a/src/operation/overlayng/OverlayNG.cpp b/src/operation/overlayng/OverlayNG.cpp
index 43f952f15..70acfbe59 100644
--- a/src/operation/overlayng/OverlayNG.cpp
+++ b/src/operation/overlayng/OverlayNG.cpp
@@ -175,7 +175,7 @@ OverlayNG::getResult()
}
else if (! inputGeom.isSingle() && inputGeom.hasPoints()) {
// handle Point-nonPoint inputs
- result = OverlayMixedPoints::overlay(opCode, ig0, ig1, pm);
+ result = OverlayMixedPoints::overlay(opCode, ig0, ig1, pm, noder);
}
else {
// handle case where both inputs are formed of edges (Lines and Polygons)
@@ -309,10 +309,10 @@ OverlayNG::extractResult(int p_opCode, OverlayGraph* graph)
//--- Build polygons
std::vector<OverlayEdge*> resultAreaEdges = graph->getResultAreaEdges();
PolygonBuilder polyBuilder(resultAreaEdges, geomFact);
- std::vector<std::unique_ptr<Polygon>> resultPolyList = polyBuilder.getPolygons();
+ std::vector<std::unique_ptr<Surface>> resultPolyList = polyBuilder.getPolygons();
bool hasResultAreaComponents = (!resultPolyList.empty());
- std::vector<std::unique_ptr<LineString>> resultLineList;
+ std::vector<std::unique_ptr<Curve>> resultLineList;
std::vector<std::unique_ptr<Point>> resultPointList;
GEOS_CHECK_FOR_INTERRUPTS();
diff --git a/src/operation/overlayng/OverlayUtil.cpp b/src/operation/overlayng/OverlayUtil.cpp
index aa7138a45..313837a34 100644
--- a/src/operation/overlayng/OverlayUtil.cpp
+++ b/src/operation/overlayng/OverlayUtil.cpp
@@ -16,6 +16,7 @@
#include <geos/operation/overlayng/InputGeometry.h>
#include <geos/operation/overlayng/OverlayGraph.h>
+#include <geos/geom/CircularString.h>
#include <geos/geom/Coordinate.h>
#include <geos/geom/Envelope.h>
#include <geos/geom/GeometryFactory.h>
@@ -279,8 +280,8 @@ OverlayUtil::isResultAreaConsistent(
/*public static*/
std::unique_ptr<Geometry>
OverlayUtil::createResultGeometry(
- std::vector<std::unique_ptr<Polygon>>& resultPolyList,
- std::vector<std::unique_ptr<LineString>>& resultLineList,
+ std::vector<std::unique_ptr<Surface>>& resultPolyList,
+ std::vector<std::unique_ptr<Curve>>& resultLineList,
std::vector<std::unique_ptr<Point>>& resultPointList,
const GeometryFactory* geometryFactory)
{
@@ -305,16 +306,22 @@ OverlayUtil::createResultGeometry(
std::unique_ptr<Geometry>
OverlayUtil::toLines(OverlayGraph* graph, bool isOutputEdges, const GeometryFactory* geomFact)
{
- std::vector<std::unique_ptr<LineString>> lines;
+ std::vector<std::unique_ptr<Geometry>> lines;
std::vector<OverlayEdge*>& edges = graph->getEdges();
for (OverlayEdge* edge : edges) {
bool includeEdge = isOutputEdges || edge->isInResultArea();
if (! includeEdge) continue;
const auto& pts = edge->getCoordinatesOriented();
- std::unique_ptr<LineString> line = geomFact->createLineString(pts);
+
+ std::unique_ptr<Geometry> edgeGeom;
+ if (edge->isCurved()) {
+ edgeGeom = geomFact->createCircularString(pts);
+ } else {
+ edgeGeom = geomFact->createLineString(pts);
+ }
// line->setUserData(labelForResult(edge));
- lines.push_back(std::move(line));
+ lines.push_back(std::move(edgeGeom));
}
return geomFact->buildGeometry(std::move(lines));
}
diff --git a/src/operation/overlayng/PolygonBuilder.cpp b/src/operation/overlayng/PolygonBuilder.cpp
index 8cc88d306..d134832a4 100644
--- a/src/operation/overlayng/PolygonBuilder.cpp
+++ b/src/operation/overlayng/PolygonBuilder.cpp
@@ -20,7 +20,7 @@
#include <geos/util/TopologyException.h>
-using geos::geom::Polygon;
+using geos::geom::Surface;
namespace geos { // geos
namespace operation { // geos.operation
@@ -28,7 +28,7 @@ namespace overlayng { // geos.operation.overlayng
/*public*/
-std::vector<std::unique_ptr<Polygon>>
+std::vector<std::unique_ptr<Surface>>
PolygonBuilder::getPolygons() const
{
return computePolygons(shellList);
@@ -42,14 +42,14 @@ PolygonBuilder::getShellRings() const
}
/*private*/
-std::vector<std::unique_ptr<Polygon>>
+std::vector<std::unique_ptr<Surface>>
PolygonBuilder::computePolygons(const std::vector<OverlayEdgeRing*>& shells) const
{
- std::vector<std::unique_ptr<Polygon>> resultPolyList;
+ std::vector<std::unique_ptr<Surface>> resultPolyList;
resultPolyList.reserve(shells.size());
// add Polygons for all shells
for (OverlayEdgeRing* er : shells) {
- std::unique_ptr<Polygon> poly = er->toPolygon(geometryFactory);
+ std::unique_ptr<Surface> poly = er->toSurface(geometryFactory);
resultPolyList.push_back(std::move(poly));
}
return resultPolyList;
diff --git a/tests/unit/operation/overlayng/OverlayNGTest.cpp b/tests/unit/operation/overlayng/OverlayNGTest.cpp
index ce7323121..92f1d32a1 100644
--- a/tests/unit/operation/overlayng/OverlayNGTest.cpp
+++ b/tests/unit/operation/overlayng/OverlayNGTest.cpp
@@ -6,12 +6,16 @@
// geos
#include <geos/operation/overlayng/OverlayNG.h>
+#include <geos/noding/SimpleNoder.h>
+#include <geos/algorithm/CircularArcIntersector.h>
+#include <geos/noding/ArcIntersectionAdder.h>
// std
#include <memory>
using namespace geos::geom;
using namespace geos::operation::overlayng;
+using namespace geos::noding;
using geos::io::WKTReader;
using geos::io::WKTWriter;
@@ -25,6 +29,13 @@ struct test_overlayng_data {
WKTReader r;
WKTWriter w;
+ double tol{0.0};
+
+ void
+ setEqualityTolerance(double d)
+ {
+ tol = d;
+ }
void
testOverlay(const std::string& a, const std::string& b, const std::string& expected, int opCode, double scaleFactor)
@@ -38,10 +49,25 @@ struct test_overlayng_data {
std::unique_ptr<Geometry> geom_a = r.read(a);
std::unique_ptr<Geometry> geom_b = r.read(b);
std::unique_ptr<Geometry> geom_expected = r.read(expected);
- std::unique_ptr<Geometry> geom_result = OverlayNG::overlay(geom_a.get(), geom_b.get(), opCode, pm.get());
- // std::string wkt_result = w.write(geom_result.get());
- // std::cout << std::endl << wkt_result << std::endl;
- ensure_equals_geometry(geom_expected.get(), geom_result.get());
+
+
+ std::unique_ptr<Geometry> geom_result;
+ if (geom_a->hasCurvedComponents() || geom_b->hasCurvedComponents()) {
+ SimpleNoder noder;
+ geos::algorithm::CircularArcIntersector cai;
+ ArcIntersectionAdder aia(cai);
+
+ noder.setArcIntersector(aia);
+ OverlayNG ong(geom_a.get(), geom_b.get(), pm.get(), opCode);
+ ong.setNoder(&noder);
+
+ geom_result = ong.getResult();
+ } else {
+ geom_result = OverlayNG::overlay(geom_a.get(), geom_b.get(), opCode, pm.get());
+ }
+ //std::string wkt_result = w.write(geom_result.get());
+ //std::cout << std::endl << wkt_result << std::endl;
+ ensure_equals_geometry_xyzm(geom_result.get(), geom_expected.get(), tol);
}
void
@@ -59,7 +85,7 @@ struct test_overlayng_data {
std::unique_ptr<Geometry> geom_result = OverlayNG::overlay(geom_a.get(), geom_b.get(), opCode, pm.get());
// std::string wkt_result = w.write(geom_result.get());
// std::cout << std::endl << wkt_result << std::endl;
- ensure_equals_exact_geometry(geom_expected.get(), geom_result.get(), 0);
+ ensure_equals_exact_geometry(geom_expected.get(), geom_result.get(), tol);
}
void
@@ -582,4 +608,197 @@ void object::test<45> ()
testOverlay(a, b, exp, OverlayNG::INTERSECTION, 0);
}
+template<>
+template<>
+void object::test<46> ()
+{
+ set_test_name("CurvePolygon/CurvePolygon intersection -> CurvePolygon");
+
+ std::string a = "CURVEPOLYGON (COMPOUNDCURVE((10 0, 0 0, 0 10, 10 10), CIRCULARSTRING (10 10, 15 5, 10 0)))";
+ std::string b = "CURVEPOLYGON (COMPOUNDCURVE((10 10, 20 10, 20 0, 10 0), CIRCULARSTRING (10 0, 5 5, 10 10)))";
+ std::string exp = "CURVEPOLYGON (CIRCULARSTRING (10 10, 15 5, 10 0, 5 5, 10 10))";
+ testOverlay(a, b, exp, OverlayNG::INTERSECTION, 0);
+}
+
+template<>
+template<>
+void object::test<47>()
+{
+ set_test_name("CurvePolygon/CurvePolygon symdifference -> MultiSurface");
+
+ std::string a = "CURVEPOLYGON (COMPOUNDCURVE((10 0, 0 0, 0 10, 10 10), CIRCULARSTRING (10 10, 15 5, 10 0)))";
+ std::string b = "CURVEPOLYGON (COMPOUNDCURVE((10 10, 20 10, 20 0, 10 0), CIRCULARSTRING (10 0, 5 5, 10 10)))";
+ std::string exp = "MULTISURFACE (CURVEPOLYGON (COMPOUNDCURVE ((10 0, 0 0, 0 10, 10 10), CIRCULARSTRING (10 10, 5 5, 10 0))), CURVEPOLYGON (COMPOUNDCURVE (CIRCULARSTRING (10 0, 15 5, 10 10), (10 10, 20 10, 20 0, 10 0))))";
+ testOverlay(a, b, exp, OverlayNG::SYMDIFFERENCE, 0);
+}
+
+template<>
+template<>
+void object::test<48>()
+{
+ set_test_name("Polygon/CircularString intersection -> CircularString");
+
+ std::string a = "POLYGON ((10 0, 10 10, 0 10, 0 0, 10 0))";
+ std::string b = "CIRCULARSTRING (5 0, 10 5, 15 0)";
+ std::string exp = "CIRCULARSTRING (5 0, 6.464466094067262 3.5355339059327378, 10 5)";
+
+ testOverlay(a, b, exp, OverlayNG::INTERSECTION, 0);
+}
+
+template<>
+template<>
+void object::test<49>()
+{
+ set_test_name("CircularString / CircularString intersection -> MultiPoint");
+
+ std::string a = "CIRCULARSTRING (-5 0, 0 5, 5 0)";
+ std::string b = "CIRCULARSTRING (-5 5, 0 0, 5 5)";
+ std::string exp = "MULTIPOINT ((4.330127018922194 2.5), (-4.330127018922194 2.5))";
+
+ testOverlay(a, b, exp, OverlayNG::INTERSECTION, 0);
+}
+
+template<>
+template<>
+void object::test<50>()
+{
+ set_test_name("CircularString / CircularString intersection -> MultiCurve");
+ std::string a = "CIRCULARSTRING (-5 0, 0 5, 5 0)";
+ std::string b = "CIRCULARSTRING (4 3, 0 -5, -4 3)";
+ std::string exp = "MULTICURVE (CIRCULARSTRING (-5 0, -4.743416490252569 1.58113883008419, -4 3), CIRCULARSTRING (4 3, 4.743416490252569 1.5811388300841898, 5 0))";
+
+ testOverlay(a, b, exp, OverlayNG::INTERSECTION, 0);
+}
+
+template<>
+template<>
+void object::test<51>()
+{
+ set_test_name("CircularString / CircularString difference -> CircularString");
+ std::string a = "CIRCULARSTRING (-5 0, 0 5, 5 0)";
+ std::string b = "CIRCULARSTRING (4 3, 0 -5, -4 3)";
+ std::string exp = "CIRCULARSTRING (-4 3, 0 5, 4 3)";
+
+ setEqualityTolerance(1e-15);
+ testOverlay(a, b, exp, OverlayNG::DIFFERENCE, 0);
+}
+
+template<>
+template<>
+void object::test<52>()
+{
+ set_test_name("CircularString / CircularString symdifference -> MultiCurve");
+ std::string a = "CIRCULARSTRING (-5 0, 0 5, 5 0)";
+ std::string b = "CIRCULARSTRING (4 3, 0 -5, -4 3)";
+ std::string exp = "MULTICURVE (CIRCULARSTRING (-4 3, 0 5, 4 3), CIRCULARSTRING (5 0, 0 -5, -5 0))";
+
+ setEqualityTolerance(1e-15);
+ testOverlay(a, b, exp, OverlayNG::SYMDIFFERENCE, 0);
+}
+
+template<>
+template<>
+void object::test<53>()
+{
+ set_test_name("CompoundCurves / LineString -> MultiPoint");
+
+ std::string a = "COMPOUNDCURVE ((-5 0, -1 0), CIRCULARSTRING (-1 0, 0 1, 1 0))";
+ std::string b = "LINESTRING (-3 -1, -4 1, 1 0)";
+ std::string exp = "MULTIPOINT ((-0.923076923076923 0.3846153846153846), (1 0), (-3.5 0))";
+
+ setEqualityTolerance(1e-15);
+ testOverlay(a, b, exp, OverlayNG::INTERSECTION, 0);
+}
+
+template<>
+template<>
+void object::test<54>()
+{
+ set_test_name("CircularString / CircularString -> Point (tangent)");
+
+ std::string a = "CIRCULARSTRING (5 3, 5 0, 0 -5, -5 0, 4 3)";
+ std::string b = "CIRCULARSTRING (-5 10, 0 15, 5 10, 4 7, -5 10)";
+ std::string exp = "POINT (0 5)";
+
+ testOverlay(a, b, exp, OverlayNG::INTERSECTION, 0);
+}
+
+template<>
+template<>
+void object::test<55>()
+{
+ set_test_name("Polygon / CurvePolygon -> CurvePolygon");
+
+ std::string a = "POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))";
+ std::string b = "CURVEPOLYGON (CIRCULARSTRING (4 5, 5 6, 6 5, 5 4, 4 5))";
+ std::string exp = "CURVEPOLYGON ((0 0, 10 0, 10 10, 0 10, 0 0), CIRCULARSTRING (4 5, 5 6, 6 5, 5 4, 4 5))";
+
+ testOverlay(a, b, exp, OverlayNG::DIFFERENCE, 0);
+}
+
+template<>
+template<>
+void object::test<56>()
+{
+ set_test_name("CircularString / Point -> Point");
+
+ std::string a = "CIRCULARSTRING (-5 0, 0 5, 5 0)";
+ std::string b = "POINT (4 3)";
+ std::string exp = b;
+
+ testOverlay(a, b, exp, OverlayNG::INTERSECTION, 0);
+}
+
+template<>
+template<>
+void object::test<57>()
+{
+ set_test_name("CompoundCurve / MultiPoint -> MultiPoint");
+
+ std::string a = "COMPOUNDCURVE ((-10 5, -5 0), CIRCULARSTRING (-5 0, 0 5, 5 0))";
+ std::string b = "MULTIPOINT (4 3, -7 2, 0 0)";
+ std::string exp = "MULTIPOINT (4 3, -7 2)";
+
+ testOverlay(a, b, exp, OverlayNG::INTERSECTION, 0);
+}
+
+template<>
+template<>
+void object::test<58>()
+{
+ set_test_name("CurvePolygon / MultiPoint -> MultiPoint");
+
+ std::string a = "CURVEPOLYGON (COMPOUNDCURVE( CIRCULARSTRING(-5 0, 0 5, 5 0), (5 0, -5 0)))";
+ std::string b = "MULTIPOINT (0 0, 4 3, 2 2, 4 4)";
+ std::string exp = "MULTIPOINT (0 0, 4 3, 2 2)";
+
+ testOverlay(a, b, exp, OverlayNG::INTERSECTION, 0);
+}
+
+template<>
+template<>
+void object::test<59>()
+{
+ set_test_name("MultiSurface / MultiPoint -> MultiPoint");
+
+ std::string a = "MULTISURFACE(CURVEPOLYGON (COMPOUNDCURVE( CIRCULARSTRING(-5 0, 0 5, 5 0), (5 0, -5 0))), ((4 4, 5 4, 5 5, 4 4)))";
+ std::string b = "MULTIPOINT (0 0, 4 3, 2 2, 4.5 4.1)";
+ std::string exp = b;
+
+ testOverlay(a, b, exp, OverlayNG::INTERSECTION, 0);
+}
+
+template<>
+template<>
+void object::test<60>()
+{
+ set_test_name("MultiSurface / MultiPoint -> MultiSurface");
+
+ std::string a = "MULTISURFACE(CURVEPOLYGON (COMPOUNDCURVE( CIRCULARSTRING(-5 0, 0 5, 5 0), (5 0, -5 0))), ((4 4, 5 4, 5 5, 4 4)))";
+ std::string b = "MULTIPOINT (0 0, 4 3, 2 2, 4.5 4.1)";
+ std::string exp = a;
+
+ testOverlay(a, b, exp, OverlayNG::UNION, 0);
+}
+
} // namespace tut
commit fa14b9a717e2231679859d20adae502ebff07492
Author: Daniel Baston <dbaston at gmail.com>
Date: Thu Mar 19 13:04:00 2026 -0400
EdgeNodingBuilder: Accept CircularString input
diff --git a/include/geos/operation/overlayng/Edge.h b/include/geos/operation/overlayng/Edge.h
index 2af2a0c67..7ec26eb16 100644
--- a/include/geos/operation/overlayng/Edge.h
+++ b/include/geos/operation/overlayng/Edge.h
@@ -64,6 +64,7 @@ private:
int bDepthDelta = 0;
bool bIsHole = false;
std::shared_ptr<const geom::CoordinateSequence> pts;
+ bool ptsCurved;
// Methods
@@ -183,6 +184,7 @@ private:
public:
+#if 0
Edge()
: aDim(OverlayLabel::DIM_UNKNOWN)
, aDepthDelta(0)
@@ -192,15 +194,20 @@ public:
, bIsHole(false)
, pts(nullptr)
{};
+#endif
+
+ Edge(const std::shared_ptr<const geom::CoordinateSequence>& p_pts, const EdgeSourceInfo* info, bool isCurved);
friend std::ostream& operator<<(std::ostream& os, const Edge& e);
static bool isCollapsed(const geom::CoordinateSequence* pts);
- Edge(const std::shared_ptr<const geom::CoordinateSequence>& p_pts, const EdgeSourceInfo* info);
+ bool isCurved() const {
+ return ptsCurved;
+ }
// return a clone of the underlying points
- std::unique_ptr<geom::CoordinateSequence> getCoordinates()
+ std::unique_ptr<geom::CoordinateSequence> getCoordinates() const
{
return pts->clone();
};
diff --git a/include/geos/operation/overlayng/EdgeNodingBuilder.h b/include/geos/operation/overlayng/EdgeNodingBuilder.h
index 0927241ba..d981a51d5 100644
--- a/include/geos/operation/overlayng/EdgeNodingBuilder.h
+++ b/include/geos/operation/overlayng/EdgeNodingBuilder.h
@@ -20,12 +20,14 @@
#include <geos/algorithm/LineIntersector.h>
#include <geos/geom/CoordinateSequence.h>
+#include <geos/geom/CircularArc.h>
#include <geos/geom/Envelope.h>
#include <geos/geom/Geometry.h>
#include <geos/geom/GeometryCollection.h>
#include <geos/geom/LinearRing.h>
#include <geos/geom/LineString.h>
#include <geos/geom/Polygon.h>
+#include <geos/geom/SimpleCurve.h>
#include <geos/noding/IntersectionAdder.h>
#include <geos/noding/Noder.h>
#include <geos/noding/SegmentString.h>
@@ -65,11 +67,16 @@ class GEOS_DLL EdgeNodingBuilder {
using PrecisionModel = geos::geom::PrecisionModel;
using Envelope = geos::geom::Envelope;
using GeometryCollection = geos::geom::GeometryCollection;
+ using Surface = geos::geom::Surface;
using Polygon = geos::geom::Polygon;
using CoordinateSequence = geos::geom::CoordinateSequence;
+ using CompoundCurve = geos::geom::CompoundCurve;
+ using CircularString = geos::geom::CircularString;
using LinearRing = geos::geom::LinearRing;
using LineString = geos::geom::LineString;
using Geometry = geos::geom::Geometry;
+ using SimpleCurve = geos::geom::SimpleCurve;
+ using Curve = geos::geom::Curve;
private:
@@ -79,7 +86,7 @@ private:
// Members
const PrecisionModel* pm;
- std::unique_ptr<std::vector<noding::SegmentString*>> inputEdges;
+ std::vector<noding::PathString*> inputEdges;
noding::Noder* customNoder;
std::array<bool, 2> hasEdges;
const Envelope* clipEnv;
@@ -112,11 +119,17 @@ private:
void addCollection(const GeometryCollection* gc, uint8_t geomIndex);
void addGeometryCollection(const GeometryCollection* gc, uint8_t geomIndex, int expectedDim);
- void addPolygon(const Polygon* poly, uint8_t geomIndex);
- void addPolygonRing(const LinearRing* ring, bool isHole, uint8_t geomIndex);
- void addLine(const LineString* line, uint8_t geomIndex);
- void addLine(std::unique_ptr<CoordinateSequence>& pts, uint8_t geomIndex);
- void addEdge(std::unique_ptr<CoordinateSequence>& cas, const EdgeSourceInfo* info);
+ void addPolygon(const Surface* poly, uint8_t geomIndex);
+ void addPolygonRing(const Curve* ring, bool isHole, uint8_t geomIndex);
+ void addSimpleCurve(const SimpleCurve* line, uint8_t geomIndex);
+ void addCompoundCurve(const CompoundCurve* line, uint8_t geomIndex);
+
+ void addCurve(const std::shared_ptr<const CoordinateSequence>& pts, const std::vector<geom::CircularArc>& arcs, uint8_t geomIndex);
+ void addLine(const std::shared_ptr<const CoordinateSequence>& pts, uint8_t geomIndex);
+
+ void addEdge(const std::shared_ptr<const CoordinateSequence>& cas, const EdgeSourceInfo* info);
+ void addCurvedEdge(const std::shared_ptr<const CoordinateSequence>& cas, const std::vector<geom::CircularArc>& arcs, const EdgeSourceInfo* info);
+
// Create a EdgeSourceInfo* owned by EdgeNodingBuilder
const EdgeSourceInfo* createEdgeSourceInfo(uint8_t index, int depthDelta, bool isHole);
@@ -156,7 +169,7 @@ private:
* @param ring the line to clip
* @return the points in the clipped line
*/
- std::unique_ptr<CoordinateSequence> clip(const LinearRing* line);
+ std::shared_ptr<const CoordinateSequence> clip(const LineString* ring) const;
/**
* Removes any repeated points from a linear component.
@@ -165,9 +178,9 @@ private:
* @param line the line to process
* @return the points of the line with repeated points removed
*/
- static std::unique_ptr<CoordinateSequence> removeRepeatedPoints(const LineString* line);
+ static std::shared_ptr<const CoordinateSequence> removeRepeatedPoints(const LineString* line);
- static int computeDepthDelta(const LinearRing* ring, bool isHole);
+ static int computeDepthDelta(const CoordinateSequence* ringCoords, bool isHole);
void add(const Geometry* g, uint8_t geomIndex);
@@ -177,8 +190,8 @@ private:
* which is used to provide source topology info to the constructed Edges
* (and is then discarded).
*/
- std::vector<Edge*> node(const std::vector<noding::SegmentString*>& segStrings);
- std::vector<Edge*> createEdges(std::vector<std::unique_ptr<noding::SegmentString>> &segStrings);
+ std::vector<Edge*> node(const std::vector<noding::PathString*>& segStrings);
+ std::vector<Edge*> createEdges(std::vector<std::unique_ptr<noding::PathString>> &segStrings);
public:
@@ -190,7 +203,6 @@ public:
*/
EdgeNodingBuilder(const PrecisionModel* p_pm, noding::Noder* p_customNoder)
: pm(p_pm)
- , inputEdges(new std::vector<noding::SegmentString*>)
, customNoder(p_customNoder)
, hasEdges{{false,false}}
, clipEnv(nullptr)
diff --git a/src/operation/overlayng/Edge.cpp b/src/operation/overlayng/Edge.cpp
index 7f231d722..2c4e8634a 100644
--- a/src/operation/overlayng/Edge.cpp
+++ b/src/operation/overlayng/Edge.cpp
@@ -27,7 +27,7 @@ using namespace geos::geom;
using geos::util::GEOSException;
/*public*/
-Edge::Edge(const std::shared_ptr<const CoordinateSequence>& p_pts, const EdgeSourceInfo* info)
+Edge::Edge(const std::shared_ptr<const CoordinateSequence>& p_pts, const EdgeSourceInfo* info, bool p_isCurved)
: aDim(OverlayLabel::DIM_UNKNOWN)
, aDepthDelta(0)
, aIsHole(false)
@@ -35,6 +35,7 @@ Edge::Edge(const std::shared_ptr<const CoordinateSequence>& p_pts, const EdgeSou
, bDepthDelta(0)
, bIsHole(false)
, pts(p_pts)
+ , ptsCurved(p_isCurved)
{
copyInfo(info);
}
diff --git a/src/operation/overlayng/EdgeNodingBuilder.cpp b/src/operation/overlayng/EdgeNodingBuilder.cpp
index b063770f1..13c123315 100644
--- a/src/operation/overlayng/EdgeNodingBuilder.cpp
+++ b/src/operation/overlayng/EdgeNodingBuilder.cpp
@@ -16,10 +16,13 @@
*
**********************************************************************/
+#include <geos/algorithm/Orientation.h>
+#include <geos/geom/CompoundCurve.h>
+#include <geos/geom/CircularString.h>
#include <geos/noding/ValidatingNoder.h>
#include <geos/noding/NodedSegmentString.h>
#include <geos/noding/MCIndexNoder.h>
-#include <geos/algorithm/Orientation.h>
+#include <geos/noding/NodableArcString.h>
#include <geos/noding/snapround/SnapRoundingNoder.h>
#include <geos/operation/overlayng/EdgeNodingBuilder.h>
#include <geos/operation/overlayng/EdgeMerger.h>
@@ -36,15 +39,17 @@ using namespace geos::geom;
using geos::operation::valid::RepeatedPointRemover;
using geos::noding::snapround::SnapRoundingNoder;
using geos::noding::Noder;
+using geos::noding::NodableArcString;
using geos::noding::MCIndexNoder;
using geos::noding::ValidatingNoder;
using geos::noding::SegmentString;
using geos::noding::NodedSegmentString;
+using geos::noding::PathString;
EdgeNodingBuilder::~EdgeNodingBuilder()
{
- for (SegmentString* ss: *inputEdges) {
+ for (auto* ss: inputEdges) {
delete ss;
}
}
@@ -108,7 +113,7 @@ EdgeNodingBuilder::build(const Geometry* geom0, const Geometry* geom1)
add(geom0, 0);
add(geom1, 1);
- std::vector<Edge*> nodedEdges = node(*inputEdges);
+ std::vector<Edge*> nodedEdges = node(inputEdges);
/**
* Merge the noded edges to eliminate duplicates.
@@ -119,14 +124,14 @@ EdgeNodingBuilder::build(const Geometry* geom0, const Geometry* geom1)
/*private*/
std::vector<Edge*>
-EdgeNodingBuilder::node(const std::vector<SegmentString*>& segStrings)
+EdgeNodingBuilder::node(const std::vector<PathString*>& segStrings)
{
std::vector<Edge*> nodedEdges;
Noder* noder = getNoder();
- noder->computeNodes(segStrings);
+ noder->computePathNodes(segStrings);
- auto nodedSS = noder->getNodedSubstrings();
+ auto nodedSS = noder->getNodedPaths();
nodedEdges = createEdges(nodedSS);
@@ -135,13 +140,15 @@ EdgeNodingBuilder::node(const std::vector<SegmentString*>& segStrings)
/*private*/
std::vector<Edge*>
-EdgeNodingBuilder::createEdges(std::vector<std::unique_ptr<SegmentString>>& segStrings)
+EdgeNodingBuilder::createEdges(std::vector<std::unique_ptr<PathString>>& segStrings)
{
std::vector<Edge*> createdEdges;
for (auto& ss : segStrings) {
const auto& pts = ss->getCoordinates();
+ const bool isCurved = dynamic_cast<const noding::ArcString*>(ss.get());
+
// don't create edges from collapsed lines
if (Edge::isCollapsed(pts.get())) continue;
@@ -150,7 +157,7 @@ EdgeNodingBuilder::createEdges(std::vector<std::unique_ptr<SegmentString>>& segS
// Record that a non-collapsed edge exists for the parent geometry
hasEdges[info->getIndex()] = true;
// Allocate the new Edge locally in a std::deque
- edgeQue.emplace_back(ss->getCoordinates(), info);
+ edgeQue.emplace_back(ss->getCoordinates(), info, isCurved);
Edge* newEdge = &(edgeQue.back());
createdEdges.push_back(newEdge);
}
@@ -179,12 +186,18 @@ EdgeNodingBuilder::add(const Geometry* g, uint8_t geomIndex)
switch (g->getGeometryTypeId())
{
case GEOS_POLYGON:
- return addPolygon(static_cast<const Polygon*>(g), geomIndex);
+ case GEOS_CURVEPOLYGON:
+ return addPolygon(static_cast<const Surface*>(g), geomIndex);
case GEOS_LINESTRING:
case GEOS_LINEARRING:
- return addLine(static_cast<const LineString*>(g), geomIndex);
+ case GEOS_CIRCULARSTRING:
+ return addSimpleCurve(static_cast<const SimpleCurve*>(g), geomIndex);
+ case GEOS_COMPOUNDCURVE:
+ return addCompoundCurve(static_cast<const CompoundCurve*>(g), geomIndex);
case GEOS_MULTILINESTRING:
case GEOS_MULTIPOLYGON:
+ case GEOS_MULTICURVE:
+ case GEOS_MULTISURFACE:
return addCollection(static_cast<const GeometryCollection*>(g), geomIndex);
case GEOS_GEOMETRYCOLLECTION:
return addGeometryCollection(static_cast<const GeometryCollection*>(g), geomIndex, g->getDimension());
@@ -192,6 +205,7 @@ EdgeNodingBuilder::add(const Geometry* g, uint8_t geomIndex)
case GEOS_MULTIPOINT:
return; // do nothing
default:
+ throw util::IllegalArgumentException("Unexpected geometry type: " + g->getGeometryType());
return; // do nothing
}
}
@@ -221,13 +235,13 @@ EdgeNodingBuilder::addGeometryCollection(const GeometryCollection* gc, uint8_t g
/*private*/
void
-EdgeNodingBuilder::addPolygon(const Polygon* poly, uint8_t geomIndex)
+EdgeNodingBuilder::addPolygon(const Surface* poly, uint8_t geomIndex)
{
- const LinearRing* shell = poly->getExteriorRing();
+ const Curve* shell = poly->getExteriorRing();
addPolygonRing(shell, false, geomIndex);
for (std::size_t i = 0; i < poly->getNumInteriorRing(); i++) {
- const LinearRing* hole = poly->getInteriorRingN(i);
+ const Curve* hole = poly->getInteriorRingN(i);
// Holes are topologically labelled opposite to the shell, since
// the interior of the polygon lies on their opposite side
@@ -238,25 +252,60 @@ EdgeNodingBuilder::addPolygon(const Polygon* poly, uint8_t geomIndex)
/*private*/
void
-EdgeNodingBuilder::addPolygonRing(const LinearRing* ring, bool isHole, uint8_t geomIndex)
- {
+EdgeNodingBuilder::addPolygonRing(const Curve* ring, bool isHole, uint8_t geomIndex)
+{
// don't add empty rings
if (ring->isEmpty()) return;
if (isClippedCompletely(ring->getEnvelopeInternal()))
return;
- std::unique_ptr<geom::CoordinateSequence> pts = clip(ring);
+ const auto type = ring->getGeometryTypeId();
+
+ if (type == GEOS_COMPOUNDCURVE) {
+ const CompoundCurve* cc = static_cast<const CompoundCurve*>(ring);
+ const auto coords = ring->getCoordinates();
+ const int depthDelta = computeDepthDelta(coords.get(), isHole);
+
+ const EdgeSourceInfo* eso = createEdgeSourceInfo(geomIndex, depthDelta, isHole);
+
+ for (std::size_t i = 0; i < cc->getNumCurves(); i++) {
+ const SimpleCurve* section = cc->getCurveN(i);
+ if (section->getNumPoints() >= 2) {
+ if (section->getGeometryTypeId() == GEOS_CIRCULARSTRING) {
+ const CircularString* cs = static_cast<const CircularString*>(section);
+ addCurvedEdge(cs->getSharedCoordinates(), cs->getArcs(), eso);
+ } else {
+ addEdge(section->getSharedCoordinates(), eso);
+ }
+ }
+ }
- /**
- * Don't add edges that collapse to a point
- */
- if (pts->size() < 2) {
return;
}
- int depthDelta = computeDepthDelta(ring, isHole);
- addEdge(pts, createEdgeSourceInfo(geomIndex, depthDelta, isHole));
+ const SimpleCurve* scRing = static_cast<const SimpleCurve*>(ring);
+ const int depthDelta = computeDepthDelta(scRing->getCoordinatesRO(), isHole);
+ const EdgeSourceInfo* esInfo = createEdgeSourceInfo(geomIndex, depthDelta, isHole);
+
+ if (ring->getGeometryTypeId() == GEOS_LINEARRING || ring->getGeometryTypeId() == GEOS_LINESTRING) {
+ // TODO: Support CircularString in RingClipper
+ std::shared_ptr<const CoordinateSequence> pts = clip(static_cast<const LineString*>(ring));
+
+ /**
+ * Don't add edges that collapse to a point
+ */
+ if (pts->size() < 2) {
+ return;
+ }
+
+ addEdge(pts, esInfo);
+ } else {
+ assert(ring->getGeometryTypeId() == GEOS_CIRCULARSTRING);
+
+ const CircularString* cs = static_cast<const CircularString*>(ring);
+ addCurvedEdge(cs->getSharedCoordinates(), cs->getArcs(), esInfo);
+ }
}
/*private*/
@@ -281,13 +330,20 @@ EdgeNodingBuilder::createEdgeSourceInfo(uint8_t index, int depthDelta, bool isHo
/*private*/
void
-EdgeNodingBuilder::addEdge(std::unique_ptr<CoordinateSequence>& cas, const EdgeSourceInfo* info)
+EdgeNodingBuilder::addEdge(const std::shared_ptr<const CoordinateSequence>& cas, const EdgeSourceInfo* info)
{
// TODO: manage these internally to EdgeNodingBuilder in a std::deque,
// since they do not have a life span longer than the EdgeNodingBuilder
// in OverlayNG::buildGraph()
- NodedSegmentString* ss = new NodedSegmentString(std::move(cas), inputHasZ, inputHasM, reinterpret_cast<const void*>(info));
- inputEdges->push_back(ss);
+ NodedSegmentString* ss = new NodedSegmentString(cas, inputHasZ, inputHasM, reinterpret_cast<const void*>(info));
+ inputEdges.push_back(ss);
+}
+
+void
+EdgeNodingBuilder::addCurvedEdge(const std::shared_ptr<const CoordinateSequence>& cas, const std::vector<CircularArc>& arcs, const EdgeSourceInfo* info)
+{
+ NodableArcString* as = new NodableArcString(arcs, cas, inputHasZ, inputHasM, reinterpret_cast<const void*>(info));
+ inputEdges.push_back(as);
}
/*private*/
@@ -299,8 +355,8 @@ EdgeNodingBuilder::isClippedCompletely(const Envelope* env) const
}
/* private */
-std::unique_ptr<geom::CoordinateSequence>
-EdgeNodingBuilder::clip(const LinearRing* ring)
+std::shared_ptr<const CoordinateSequence>
+EdgeNodingBuilder::clip(const LineString* ring) const
{
const Envelope* env = ring->getEnvelopeInternal();
@@ -316,16 +372,20 @@ EdgeNodingBuilder::clip(const LinearRing* ring)
}
/*private*/
-std::unique_ptr<CoordinateSequence>
+std::shared_ptr<const CoordinateSequence>
EdgeNodingBuilder::removeRepeatedPoints(const LineString* line)
{
- const CoordinateSequence* pts = line->getCoordinatesRO();
- return RepeatedPointRemover::removeRepeatedPoints(pts);
+ const std::shared_ptr<const CoordinateSequence>& pts = line->getSharedCoordinates();
+ if (pts->hasRepeatedPoints()) {
+ return RepeatedPointRemover::removeRepeatedPoints(pts.get());
+ } else {
+ return pts;
+ }
}
/*private*/
int
-EdgeNodingBuilder::computeDepthDelta(const LinearRing* ring, bool isHole)
+EdgeNodingBuilder::computeDepthDelta(const CoordinateSequence* ringCoords, bool isHole)
{
/**
* Compute the orientation of the ring, to
@@ -335,7 +395,8 @@ EdgeNodingBuilder::computeDepthDelta(const LinearRing* ring, bool isHole)
* It is important to compute orientation on the original ring,
* since topology collapse can make the orientation computation give the wrong answer.
*/
- bool isCCW = algorithm::Orientation::isCCW(ring->getCoordinatesRO());
+
+ const bool isCCW = algorithm::Orientation::isCCW(ringCoords);
/**
* Compute whether ring is in canonical orientation or not.
@@ -360,29 +421,46 @@ EdgeNodingBuilder::computeDepthDelta(const LinearRing* ring, bool isHole)
/*private*/
void
-EdgeNodingBuilder::addLine(const LineString* line, uint8_t geomIndex)
+EdgeNodingBuilder::addSimpleCurve(const SimpleCurve* curve, uint8_t geomIndex)
{
// don't add empty lines
- if (line->isEmpty()) return;
+ if (curve->isEmpty()) return;
- if (isClippedCompletely(line->getEnvelopeInternal()))
+ if (isClippedCompletely(curve->getEnvelopeInternal()))
return;
+ if (curve->getGeometryTypeId() == GEOS_CIRCULARSTRING) {
+ // TODO: Support CircularString in LineLimiter
+ const CircularString* cs = static_cast<const CircularString*>(curve);
+ addCurve(cs->getSharedCoordinates(), cs->getArcs(), geomIndex);
+ return;
+ }
+
+ const LineString* line = detail::down_cast<const LineString*>(curve);
+
if (isToBeLimited(line)) {
std::vector<std::unique_ptr<CoordinateSequence>>& sections = limit(line);
for (auto& pts : sections) {
- addLine(pts, geomIndex);
+ addLine(std::move(pts), geomIndex);
}
}
else {
- std::unique_ptr<CoordinateSequence> ptsNoRepeat = removeRepeatedPoints(line);
+ const auto ptsNoRepeat = removeRepeatedPoints(line);
addLine(ptsNoRepeat, geomIndex);
}
}
+void
+EdgeNodingBuilder::addCompoundCurve(const CompoundCurve* curve, uint8_t geomIndex)
+{
+ for (std::size_t i = 0; i < curve->getNumCurves(); i++) {
+ addSimpleCurve(curve->getCurveN(i), geomIndex);
+ }
+}
+
/*private*/
void
-EdgeNodingBuilder::addLine(std::unique_ptr<CoordinateSequence>& pts, uint8_t geomIndex)
+EdgeNodingBuilder::addLine(const std::shared_ptr<const CoordinateSequence>& pts, uint8_t geomIndex)
{
/**
* Don't add edges that collapse to a point
@@ -394,6 +472,20 @@ EdgeNodingBuilder::addLine(std::unique_ptr<CoordinateSequence>& pts, uint8_t geo
addEdge(pts, createEdgeSourceInfo(geomIndex));
}
+/*private*/
+void
+EdgeNodingBuilder::addCurve(const std::shared_ptr<const CoordinateSequence>& pts, const std::vector<CircularArc>& arcs, uint8_t geomIndex)
+{
+ /**
+ * Don't add edges that collapse to a point
+ */
+ if (pts->size() < 2) {
+ return;
+ }
+
+ addCurvedEdge(pts, arcs, createEdgeSourceInfo(geomIndex));
+}
+
/*private*/
bool
EdgeNodingBuilder::isToBeLimited(const LineString* line) const
diff --git a/tests/unit/operation/overlayng/EdgeNodingBuilderTest.cpp b/tests/unit/operation/overlayng/EdgeNodingBuilderTest.cpp
new file mode 100644
index 000000000..017dffc9a
--- /dev/null
+++ b/tests/unit/operation/overlayng/EdgeNodingBuilderTest.cpp
@@ -0,0 +1,110 @@
+//
+// Test Suite for geos::operation::overlayng::EdgeNodingBuilder class.
+
+#include <tut/tut.hpp>
+#include <utility.h>
+
+// geos
+#include <geos/geom/CircularString.h>
+#include <geos/algorithm/CircularArcIntersector.h>
+#include <geos/noding/ArcIntersectionAdder.h>
+#include <geos/noding/SimpleNoder.h>
+#include <geos/operation/overlayng/EdgeNodingBuilder.h>
+
+// std
+#include <memory>
+
+
+
+using namespace geos::geom;
+using namespace geos::noding;
+using namespace geos::operation::overlayng;
+using geos::algorithm::CircularArcIntersector;
+using geos::io::WKTReader;
+using geos::io::WKTWriter;
+
+namespace tut {
+struct test_edgenodingbuilder_data {
+ WKTReader reader_;
+
+ PrecisionModel pm_floating{PrecisionModel::Type::FLOATING};
+ GeometryFactory::Ptr factory_{GeometryFactory::create(&pm_floating)};
+
+ std::unique_ptr<Geometry> toGeometry(const Edge& edge) const
+ {
+ if (edge.isCurved()) {
+ return factory_->createCircularString(edge.getCoordinates());
+ }
+ return factory_->createLineString(edge.getCoordinates());
+ }
+
+ std::unique_ptr<Geometry> toGeometry(const std::vector<Edge*> edges) const
+ {
+ std::vector<std::unique_ptr<Geometry>> geoms;
+ for (const auto* edge : edges) {
+ geoms.push_back(toGeometry(*edge));
+ }
+ return factory_->createGeometryCollection(std::move(geoms));
+ }
+
+ void checkEdges(const std::string& wkt1, const std::string& wkt2, const std::string& wktExpected) const
+ {
+ auto geom1 = reader_.read(wkt1);
+ auto geom2 = reader_.read(wkt2);
+
+ SimpleNoder noder;
+ CircularArcIntersector cai;
+ ArcIntersectionAdder aia(cai);
+ noder.setArcIntersector(aia);
+
+ EdgeNodingBuilder builder(geom1->getPrecisionModel(), &noder);
+
+ auto edges = builder.build(geom1.get(), geom2.get());
+ std::unique_ptr<Geometry> expected = reader_.read(wktExpected);
+ ensure_equals_geometry_xyzm(toGeometry(edges).get(), expected.get(), 1e-5);
+ }
+};
+
+typedef test_group<test_edgenodingbuilder_data> group;
+typedef group::object object;
+
+group test_edgenodingbuilder_group("geos::operation::overlayng::EdgeNodingBuilder");
+
+template<>
+template<>
+void object::test<1>()
+{
+ set_test_name("two CircularStrings");
+
+ checkEdges(
+ "CIRCULARSTRING (0 0, 1 1, 2 0)",
+ "CIRCULARSTRING (0 1, 1 0, 2 1)",
+ "GEOMETRYCOLLECTION ("
+ "CIRCULARSTRING (0 0, 0.0340742 0.258819, 0.133975 0.5),"
+ "CIRCULARSTRING (0.133975 0.5, 1 1, 1.86603 0.5),"
+ "CIRCULARSTRING (1.86603 0.5, 1.96593 0.258819, 2 0),"
+ "CIRCULARSTRING (0 1, 0.0340742 0.741181, 0.133975 0.5),"
+ "CIRCULARSTRING (0.133975 0.5, 1 0, 1.86603 0.5),"
+ "CIRCULARSTRING (1.86603 0.5, 1.96593 0.741181, 2 1)"
+ ")");
+}
+
+template<>
+template<>
+void object::test<2>()
+{
+ set_test_name("CurvePolygon and LineString");
+
+ checkEdges("CURVEPOLYGON (COMPOUNDCURVE(CIRCULARSTRING (0 0, 1 1, 2 0), (2 0, 0 0)))",
+ "LINESTRING (1 0, 2 1)",
+ "GEOMETRYCOLLECTION ("
+ "CIRCULARSTRING (0 0, 0.6173165676349103 0.9238795325112867, 1.7071067811865475 0.7071067811865475),"
+ "CIRCULARSTRING (1.7071067811865475 0.7071067811865475, 1.9238795325112867 0.3826834323650898, 2 0),"
+ "LINESTRING (1 0, 1.7071067811865475 0.7071067811865475),"
+ "LINESTRING (1.7071067811865475 0.7071067811865475, 2 1),"
+ "LINESTRING (0 0, 1 0),"
+ "LINESTRING (1 0, 2 0)"
+ ")");
+}
+
+}
diff --git a/tests/unit/operation/overlayng/OverlayGraphTest.cpp b/tests/unit/operation/overlayng/OverlayGraphTest.cpp
index e179702fe..7fce5ff72 100644
--- a/tests/unit/operation/overlayng/OverlayGraphTest.cpp
+++ b/tests/unit/operation/overlayng/OverlayGraphTest.cpp
@@ -52,7 +52,7 @@ struct test_overlaygraph_data {
std::unique_ptr<CoordinateSequence> cs = line->getCoordinates();
EdgeSourceInfo esi(0);
- Edge e(std::move(cs), &esi);
+ Edge e(std::move(cs), &esi, false);
std::unique_ptr<CoordinateSequence> pts = e.getCoordinates();
graph->addEdge(&e);
-----------------------------------------------------------------------
Summary of changes:
include/geos/algorithm/CircularArcs.h | 15 ++
include/geos/algorithm/PointLocator.h | 6 +
include/geos/geom/CircularArc.h | 2 +
include/geos/geom/Curve.h | 2 +
include/geos/operation/overlayng/Edge.h | 42 ++--
.../geos/operation/overlayng/EdgeNodingBuilder.h | 43 +++-
include/geos/operation/overlayng/LineBuilder.h | 8 +-
include/geos/operation/overlayng/OverlayEdge.h | 10 +-
include/geos/operation/overlayng/OverlayEdgeRing.h | 23 +-
include/geos/operation/overlayng/OverlayGraph.h | 4 +-
.../geos/operation/overlayng/OverlayMixedPoints.h | 17 +-
include/geos/operation/overlayng/OverlayUtil.h | 17 +-
include/geos/operation/overlayng/PolygonBuilder.h | 6 +-
src/algorithm/CircularArcIntersector.cpp | 162 +++++++-----
src/algorithm/CircularArcs.cpp | 94 ++++---
src/algorithm/PointLocator.cpp | 80 +++++-
src/geom/CircularArc.cpp | 14 ++
src/geom/Curve.cpp | 6 +
src/geom/GeometryFactory.cpp | 22 ++
src/noding/NodableArcString.cpp | 243 +++++++++++++-----
src/noding/snap/SnappingNoder.cpp | 3 +-
src/operation/overlayng/Edge.cpp | 3 +-
src/operation/overlayng/EdgeNodingBuilder.cpp | 191 +++++++++++----
.../overlayng/IndexedPointOnLineLocator.cpp | 2 +-
src/operation/overlayng/InputGeometry.cpp | 20 +-
src/operation/overlayng/LineBuilder.cpp | 33 ++-
src/operation/overlayng/OverlayEdgeRing.cpp | 150 ++++++++----
src/operation/overlayng/OverlayGraph.cpp | 24 +-
src/operation/overlayng/OverlayMixedPoints.cpp | 40 +--
src/operation/overlayng/OverlayNG.cpp | 6 +-
src/operation/overlayng/OverlayNGRobust.cpp | 3 -
src/operation/overlayng/OverlayUtil.cpp | 39 ++-
src/operation/overlayng/PolygonBuilder.cpp | 10 +-
.../unit/algorithm/CircularArcIntersectorTest.cpp | 60 +++++
tests/unit/capi/GEOSDifferenceTest.cpp | 8 +-
tests/unit/capi/GEOSIntersectionTest.cpp | 10 +-
tests/unit/capi/GEOSNodeTest.cpp | 8 +-
tests/unit/capi/GEOSSymDifferenceTest.cpp | 10 +-
tests/unit/capi/GEOSUnionTest.cpp | 10 +-
tests/unit/geom/CircularStringTest.cpp | 8 +-
tests/unit/geom/CompoundCurveTest.cpp | 9 +-
tests/unit/geom/CurvePolygonTest.cpp | 8 +-
tests/unit/geom/MultiCurveTest.cpp | 9 +-
tests/unit/geom/MultiSurfaceTest.cpp | 8 +-
.../operation/overlayng/EdgeNodingBuilderTest.cpp | 110 +++++++++
.../unit/operation/overlayng/OverlayGraphTest.cpp | 2 +-
tests/unit/operation/overlayng/OverlayNGTest.cpp | 271 ++++++++++++++++++++-
47 files changed, 1453 insertions(+), 418 deletions(-)
create mode 100644 tests/unit/operation/overlayng/EdgeNodingBuilderTest.cpp
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list