[geos-commits] [SCM] GEOS branch main-relate-ng updated. 8ee1ae84dfdfe26101b74271fd3c426c8593df56
git at osgeo.org
git at osgeo.org
Thu Jul 25 08:10:50 PDT 2024
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-relate-ng has been updated
via 8ee1ae84dfdfe26101b74271fd3c426c8593df56 (commit)
from 4faaea0f46da942c0e3a0db9f6a5f96aa0cb8e17 (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 8ee1ae84dfdfe26101b74271fd3c426c8593df56
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Thu Jul 25 08:09:54 2024 -0700
Add EdgeSegmentIntersector EdgeSegmentOverlapAction EdgeSetIntersector TopologyComputer
diff --git a/include/geos/index/chain/MonotoneChain.h b/include/geos/index/chain/MonotoneChain.h
index 4b0dbb427..d6b1b1ad6 100644
--- a/include/geos/index/chain/MonotoneChain.h
+++ b/include/geos/index/chain/MonotoneChain.h
@@ -151,6 +151,9 @@ public:
return context;
}
+ void setId(int p_id) { id = p_id; }
+ int getId() const { return id; }
+
private:
void computeSelect(const geom::Envelope& searchEnv,
@@ -197,6 +200,10 @@ private:
/// Owned by this class
mutable geom::Envelope env;
+
+ /// Useful for optimizing chain comparisons
+ int id = 0;
+
};
} // namespace geos::index::chain
diff --git a/include/geos/operation/relateng/EdgeSegmentIntersector.h b/include/geos/operation/relateng/EdgeSegmentIntersector.h
new file mode 100644
index 000000000..5153b68fd
--- /dev/null
+++ b/include/geos/operation/relateng/EdgeSegmentIntersector.h
@@ -0,0 +1,80 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (c) 2024 Martin Davis
+ * Copyright (C) 2024 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#pragma once
+
+
+#include <geos/noding/SegmentIntersector.h>
+#include <geos/algorithm/LineIntersector.h>
+#include <geos/export.h>
+
+
+// Forward declarations
+namespace geos {
+namespace noding {
+ class SegmentString;
+}
+namespace operation {
+namespace relateng {
+ class RelateSegmentString;
+ class TopologyComputer;
+}
+}
+}
+
+
+using geos::noding::SegmentIntersector;
+using geos::noding::SegmentString;
+using geos::algorithm::LineIntersector;
+
+
+namespace geos { // geos.
+namespace operation { // geos.operation
+namespace relateng { // geos.operation.relateng
+
+class GEOS_DLL EdgeSegmentIntersector : public SegmentIntersector {
+
+private:
+
+ // Members
+ LineIntersector li;
+ TopologyComputer& topoComputer;
+
+ // Methods
+
+
+ void addIntersections(
+ RelateSegmentString* ssA, std::size_t segIndexA,
+ RelateSegmentString* ssB, std::size_t segIndexB);
+
+
+public:
+
+ EdgeSegmentIntersector(TopologyComputer& p_topoBuilder)
+ : topoComputer(p_topoBuilder)
+ {};
+
+ void processIntersections(
+ SegmentString* ss0, std::size_t segIndex0,
+ SegmentString* ss1, std::size_t segIndex1) override;
+
+ bool isDone() const override;
+
+};
+
+} // namespace geos.operation.relateng
+} // namespace geos.operation
+} // namespace geos
+
diff --git a/include/geos/operation/relateng/EdgeSegmentOverlapAction.h b/include/geos/operation/relateng/EdgeSegmentOverlapAction.h
new file mode 100644
index 000000000..7f84ec797
--- /dev/null
+++ b/include/geos/operation/relateng/EdgeSegmentOverlapAction.h
@@ -0,0 +1,75 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2001-2002 Vivid Solutions Inc.
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************
+ *
+ * Last port: index/chain/MonotoneChainOverlapAction.java rev. 1.6 (JTS-1.10)
+ *
+ **********************************************************************/
+
+#pragma once
+
+#include <geos/export.h>
+#include <geos/geom/LineSegment.h>
+#include <geos/index/chain/MonotoneChainOverlapAction.h>
+
+
+// Forward declarations
+namespace geos {
+namespace index {
+namespace chain {
+ class MonotoneChain;
+}
+}
+namespace noding {
+ class SegmentIntersector;
+}
+}
+
+
+using geos::index::chain::MonotoneChain;
+using geos::index::chain::MonotoneChainOverlapAction;
+using geos::noding::SegmentIntersector;
+
+
+namespace geos {
+namespace operation { // geos::operation
+namespace relateng { // geos::operation::relateng
+
+/** \brief
+ * The action for the internal iterator for performing
+ * overlap queries on a MonotoneChain.
+ */
+class GEOS_DLL EdgeSegmentOverlapAction : public MonotoneChainOverlapAction {
+
+private:
+
+ SegmentIntersector& si;
+
+
+public:
+
+ EdgeSegmentOverlapAction(SegmentIntersector& p_si)
+ : si(p_si)
+ {}
+
+ void overlap(
+ const MonotoneChain& mc1, std::size_t start1,
+ const MonotoneChain& mc2, std::size_t start2) override;
+
+
+};
+
+} // namespace geos::index::chain
+} // namespace geos::index
+} // namespace geos
+
diff --git a/include/geos/operation/relateng/EdgeSetIntersector.h b/include/geos/operation/relateng/EdgeSetIntersector.h
new file mode 100644
index 000000000..bf95bb6d7
--- /dev/null
+++ b/include/geos/operation/relateng/EdgeSetIntersector.h
@@ -0,0 +1,95 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (c) 2024 Martin Davis
+ * Copyright (C) 2024 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#pragma once
+
+
+#include <geos/index/strtree/TemplateSTRtree.h>
+#include <geos/index/chain/MonotoneChain.h>
+#include <geos/export.h>
+
+
+// Forward declarations
+namespace geos {
+namespace geom {
+ class Geometry;
+ class Envelope;
+}
+namespace noding {
+ class SegmentString;
+}
+namespace operation {
+namespace relateng {
+ class RelateSegmentString;
+ class EdgeSegmentIntersector;
+}
+}
+}
+
+
+using geos::geom::Envelope;
+using geos::geom::Geometry;
+using geos::index::strtree::TemplateSTRtree;
+using geos::index::chain::MonotoneChain;
+using geos::operation::relateng::EdgeSegmentIntersector;
+
+
+namespace geos { // geos.
+namespace operation { // geos.operation
+namespace relateng { // geos.operation.relateng
+
+class GEOS_DLL EdgeSetIntersector {
+
+private:
+
+ // Members
+ TemplateSTRtree<const MonotoneChain*> index;
+ // HPRtree index = new HPRtree();
+ const Envelope* envelope = nullptr;
+ std::deque<MonotoneChain> monoChains;
+ int idCounter = 0;
+
+
+ // Methods
+
+ void addToIndex(std::unique_ptr<RelateSegmentString>& segStr);
+
+ void addEdges(std::vector<std::unique_ptr<RelateSegmentString>>& segStrings);
+
+
+public:
+
+ EdgeSetIntersector(
+ std::vector<std::unique_ptr<RelateSegmentString>>& edgesA,
+ std::vector<std::unique_ptr<RelateSegmentString>>& edgesB,
+ const Envelope* env)
+ : envelope(env)
+ , idCounter(0)
+ {
+ addEdges(edgesA);
+ addEdges(edgesB);
+ // build index to ensure thread-safety
+ // index.build();
+ };
+
+ void process(EdgeSegmentIntersector& intersector);
+
+
+};
+
+} // namespace geos.operation.relateng
+} // namespace geos.operation
+} // namespace geos
+
diff --git a/include/geos/operation/relateng/NodeSection.h b/include/geos/operation/relateng/NodeSection.h
index c464d7fed..8fe00870e 100644
--- a/include/geos/operation/relateng/NodeSection.h
+++ b/include/geos/operation/relateng/NodeSection.h
@@ -15,14 +15,15 @@
#pragma once
+#include <geos/export.h>
+#include <geos/geom/Coordinate.h>
+
#include <string>
#include <sstream>
-#include <geos/export.h>
// Forward declarations
namespace geos {
namespace geom {
- class CoordinateXY;
class Geometry;
}
}
@@ -63,7 +64,7 @@ private:
const Geometry* m_poly;
bool m_isNodeAtVertex;
const CoordinateXY* m_v0;
- const CoordinateXY* m_nodePt;
+ const CoordinateXY m_nodePt;
const CoordinateXY* m_v1;
// Methods
@@ -82,7 +83,7 @@ public:
const Geometry* poly,
bool isNodeAtVertex,
const CoordinateXY* v0,
- const CoordinateXY* nodePt,
+ const CoordinateXY nodePt,
const CoordinateXY* v1)
: m_isA(isA)
, m_dim(dim)
@@ -109,7 +110,7 @@ public:
const CoordinateXY* getVertex(int i) const;
- const CoordinateXY* nodePt() const;
+ const CoordinateXY& nodePt() const;
int dimension() const;
@@ -129,6 +130,8 @@ public:
bool isArea() const;
+ static bool isAreaArea(const NodeSection& a, const NodeSection& b);
+
bool isA() const;
bool isSameGeometry(const NodeSection& ns) const;
diff --git a/include/geos/operation/relateng/RelateSegmentString.h b/include/geos/operation/relateng/RelateSegmentString.h
index f6f284e40..d16ddb35b 100644
--- a/include/geos/operation/relateng/RelateSegmentString.h
+++ b/include/geos/operation/relateng/RelateSegmentString.h
@@ -144,7 +144,7 @@ public:
const Geometry* getPolygonal() const;
- NodeSection* createNodeSection(std::size_t segIndex, const CoordinateXY& intPt) const;
+ NodeSection* createNodeSection(std::size_t segIndex, const CoordinateXY intPt) const;
/**
* Tests if a segment intersection point has that segment as its
diff --git a/include/geos/operation/relateng/TopologyComputer.h b/include/geos/operation/relateng/TopologyComputer.h
new file mode 100644
index 000000000..f4709940b
--- /dev/null
+++ b/include/geos/operation/relateng/TopologyComputer.h
@@ -0,0 +1,226 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (c) 2024 Martin Davis
+ * Copyright (C) 2024 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#pragma once
+
+#include <geos/operation/relateng/NodeSections.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/Location.h>
+#include <geos/export.h>
+
+// Forward declarations
+namespace geos {
+namespace operation {
+namespace relateng {
+ class NodeSection;
+ class RelateGeometry;
+ class RelateNode;
+ class TopologyPredicate;
+}
+}
+}
+
+
+using geos::geom::CoordinateXY;
+using geos::geom::Location;
+using geos::operation::relateng::NodeSection;
+using geos::operation::relateng::NodeSections;
+using geos::operation::relateng::RelateGeometry;
+using geos::operation::relateng::RelateNode;
+using geos::operation::relateng::TopologyPredicate;
+
+
+namespace geos { // geos.
+namespace operation { // geos.operation
+namespace relateng { // geos.operation.relateng
+
+
+class GEOS_DLL TopologyComputer {
+
+private:
+
+ // Members
+ TopologyPredicate* predicate;
+ RelateGeometry* geomA;
+ RelateGeometry* geomB;
+ std::map<CoordinateXY, NodeSections*> nodeMap;
+ std::deque<std::unique_ptr<NodeSections>> nodeSectionsStore;
+
+ // Methods
+
+ /**
+ * Determine a priori partial EXTERIOR topology based on dimensions.
+ */
+ void initExteriorDims();
+
+ void initExteriorEmpty(bool geomNonEmpty);
+
+ RelateGeometry* getGeometry(bool isA) const;
+
+ void updateDim(Location locA, Location locB, int dimension);
+
+ void updateDim(bool isAB, Location loc1, Location loc2, int dimension);
+
+ /**
+ * Update topology for an intersection between A and B.
+ *
+ * @param a the section for geometry A
+ * @param b the section for geometry B
+ */
+ void updateIntersectionAB(const NodeSection* a, const NodeSection* b);
+
+ /**
+ * Updates topology for an AB Area-Area crossing node.
+ * Sections cross at a node if (a) the intersection is proper
+ * (i.e. in the interior of two segments)
+ * or (b) if non-proper then whether the linework crosses
+ * is determined by the geometry of the segments on either side of the node.
+ * In these situations the area geometry interiors intersect (in dimension 2).
+ *
+ * @param a the section for geometry A
+ * @param b the section for geometry B
+ */
+ void updateAreaAreaCross(const NodeSection* a, const NodeSection* b);
+
+ /**
+ * Updates topology for a node at an AB edge intersection.
+ *
+ * @param a the section for geometry A
+ * @param b the section for geometry B
+ */
+ void updateNodeLocation(const NodeSection* a, const NodeSection* b);
+
+ void addNodeSections(NodeSection* ns0, NodeSection* ns1);
+
+ void addLineEndOnPoint(bool isLineA, Location locLineEnd, Location locPoint);
+
+ void addLineEndOnLine(bool isLineA, Location locLineEnd, Location locLine);
+
+ void addLineEndOnArea(bool isLineA, Location locLineEnd, Location locArea);
+
+ /**
+ * Updates topology for an area vertex (in Interior or on Boundary)
+ * intersecting a point.
+ * Note that because the largest dimension of intersecting target is determined,
+ * the intersecting point is not part of any other target geometry,
+ * and hence its neighbourhood is in the Exterior of the target.
+ *
+ * @param isAreaA whether the area is the A input
+ * @param locArea the location of the vertex in the area
+ * @param pt the point at which topology is being updated
+ */
+ void addAreaVertexOnPoint(bool isAreaA, Location locArea);
+
+ void addAreaVertexOnLine(bool isAreaA, Location locArea, Location locTarget);
+
+ void evaluateNode(NodeSections* nodeSections);
+
+ void evaluateNodeEdges(const RelateNode* node);
+
+ NodeSections* getNodeSections(const CoordinateXY& nodePt);
+
+
+
+public:
+
+ TopologyComputer(
+ TopologyPredicate* p_predicate,
+ RelateGeometry* p_geomA,
+ RelateGeometry* p_geomB)
+ : predicate(p_predicate)
+ , geomA(p_geomA)
+ , geomB(p_geomB)
+ {
+ initExteriorDims();
+ };
+
+ int getDimension(bool isA) const;
+
+ bool isAreaArea() const;
+
+ /**
+ * Indicates whether the input geometries require self-noding
+ * for correct evaluation of specific spatial predicates.
+ * Self-noding is required for geometries which may self-cross
+ * - i.e. lines, and overlapping polygons in GeometryCollections.
+ * Self-noding is not required for polygonal geometries.
+ * This ensures that the locations of nodes created by
+ * crossing segments are computed explicitly.
+ * This ensures that node locations match in situations
+ * where a self-crossing and mutual crossing occur at the same logical location.
+ * E.g. a self-crossing line tested against a single segment
+ * identical to one of the crossed segments.
+ *
+ * @return true if self-noding is required
+ */
+ bool isSelfNodingRequired() const;
+
+ bool isExteriorCheckRequired(bool isA) const;
+
+ bool isResultKnown() const;
+
+ bool getResult() const;
+
+ /**
+ * Finalize the evaluation.
+ */
+ void finish();
+
+ void addIntersection(NodeSection* a, NodeSection* b);
+
+ void addPointOnPointInterior();
+
+ void addPointOnPointExterior(bool isGeomA);
+
+ void addPointOnGeometry(bool isA, Location locTarget, int dimTarget);
+
+ void addLineEndOnGeometry(bool isLineA, Location locLineEnd, Location locTarget, int dimTarget);
+
+ /**
+ * Adds topology for an area vertex interaction with a target geometry element.
+ * Assumes the target geometry element has highest dimension
+ * (i.e. if the point lies on two elements of different dimension,
+ * the location on the higher dimension element is provided.
+ * This is the semantic provided by {@link RelatePointLocator}.
+ *
+ * Note that in a GeometryCollection containing overlapping or adjacent polygons,
+ * the area vertex location may be INTERIOR instead of BOUNDARY.
+ *
+ * @param isAreaA the input that is the area
+ * @param locArea the location on the area
+ * @param locTarget the location on the target geometry element
+ * @param dimTarget the dimension of the target geometry element
+ * @param pt the point of interaction
+ */
+ void addAreaVertex(bool isAreaA, Location locArea, Location locTarget, int dimTarget);
+
+ void addAreaVertexOnArea(bool isAreaA, Location locArea, Location locTarget);
+
+ void evaluateNodes();
+
+ /**
+ * Disable copy construction and assignment. Apparently needed to make this
+ * class compile under MSVC. (See https://stackoverflow.com/q/29565299)
+ */
+ TopologyComputer(const TopologyComputer&) = delete;
+ TopologyComputer& operator=(const TopologyComputer&) = delete;
+
+
+};
+
+} // namespace geos.operation.relateng
+} // namespace geos.operation
+} // namespace geos
+
diff --git a/src/operation/relateng/AdjacentEdgeLocator.cpp b/src/operation/relateng/AdjacentEdgeLocator.cpp
index a065ee984..aeb5f8084 100644
--- a/src/operation/relateng/AdjacentEdgeLocator.cpp
+++ b/src/operation/relateng/AdjacentEdgeLocator.cpp
@@ -94,7 +94,7 @@ AdjacentEdgeLocator::createSection(const CoordinateXY* p,
if (prev->distance(*p) == 0 || next->distance(*p) == 0) {
//System.out.println("Found zero-length section segment");
};
- return new NodeSection(true, Dimension::A, 1, 0, nullptr, false, prev, p, next);
+ return new NodeSection(true, Dimension::A, 1, 0, nullptr, false, prev, *p, next);
}
diff --git a/src/operation/relateng/DimensionLocation.cpp b/src/operation/relateng/DimensionLocation.cpp
index fd5fceb4f..50a3457d8 100644
--- a/src/operation/relateng/DimensionLocation.cpp
+++ b/src/operation/relateng/DimensionLocation.cpp
@@ -37,7 +37,6 @@ DimensionLocation::locationArea(Location loc)
default:
return EXTERIOR;
}
- return EXTERIOR;
}
@@ -51,7 +50,6 @@ DimensionLocation::locationLine(Location loc)
default:
return EXTERIOR;
}
- return EXTERIOR;
}
@@ -64,7 +62,6 @@ DimensionLocation::locationPoint(Location loc)
default:
return EXTERIOR;
}
- return EXTERIOR;
}
@@ -83,7 +80,6 @@ DimensionLocation::location(int dimLoc)
default:
return Location::EXTERIOR;
}
- return Location::EXTERIOR;
}
@@ -103,7 +99,6 @@ DimensionLocation::dimension(int dimLoc)
default:
return Dimension::False;
}
- return Dimension::False;
}
@@ -113,7 +108,8 @@ DimensionLocation::dimension(int dimLoc, int exteriorDim)
{
if (dimLoc == EXTERIOR)
return exteriorDim;
- return dimension(dimLoc);
+ else
+ return dimension(dimLoc);
}
diff --git a/src/operation/relateng/EdgeSegmentIntersector.cpp b/src/operation/relateng/EdgeSegmentIntersector.cpp
new file mode 100644
index 000000000..8ea3c5205
--- /dev/null
+++ b/src/operation/relateng/EdgeSegmentIntersector.cpp
@@ -0,0 +1,114 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (c) 2024 Martin Davis
+ * Copyright (C) 2024 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#include <geos/geom/Geometry.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/operation/relateng/EdgeSegmentIntersector.h>
+#include <geos/operation/relateng/RelateSegmentString.h>
+#include <geos/operation/relateng/TopologyComputer.h>
+#include <geos/operation/relateng/NodeSection.h>
+#include <geos/index/chain/MonotoneChain.h>
+#include <geos/index/chain/MonotoneChainBuilder.h>
+
+
+using geos::geom::CoordinateXYZM;
+using geos::geom::CoordinateXY;
+using geos::geom::Geometry;
+using geos::noding::SegmentString;
+using geos::index::chain::MonotoneChain;
+using geos::index::chain::MonotoneChainBuilder;
+
+
+namespace geos { // geos
+namespace operation { // geos.operation
+namespace relateng { // geos.operation.relateng
+
+
+/* public override */
+bool
+EdgeSegmentIntersector::isDone() const
+{
+ return topoComputer.isResultKnown();
+}
+
+
+/* public override */
+void
+EdgeSegmentIntersector::processIntersections(
+ SegmentString* ss0, std::size_t segIndex0,
+ SegmentString* ss1, std::size_t segIndex1)
+{
+ // don't intersect a segment with itself
+ if (ss0 == ss1 && segIndex0 == segIndex1)
+ return;
+
+ RelateSegmentString* rss0 = static_cast<RelateSegmentString*>(ss0);
+ RelateSegmentString* rss1 = static_cast<RelateSegmentString*>(ss1);
+ //TODO: move this ordering logic to TopologyBuilder
+ if (rss0->isA()) {
+ addIntersections(rss0, segIndex0, rss1, segIndex1);
+ }
+ else {
+ addIntersections(rss1, segIndex1, rss0, segIndex0);
+ }
+}
+
+
+/* private */
+void
+EdgeSegmentIntersector::addIntersections(
+ RelateSegmentString* ssA, std::size_t segIndexA,
+ RelateSegmentString* ssB, std::size_t segIndexB)
+{
+
+ const CoordinateXY& a0 = ssA->getCoordinate(segIndexA);
+ const CoordinateXY& a1 = ssA->getCoordinate(segIndexA + 1);
+ const CoordinateXY& b0 = ssB->getCoordinate(segIndexB);
+ const CoordinateXY& b1 = ssB->getCoordinate(segIndexB + 1);
+
+ li.computeIntersection(a0, a1, b0, b1);
+
+ if (! li.hasIntersection())
+ return;
+
+ for (uint32_t i = 0; i < li.getIntersectionNum(); i++)
+ {
+ const CoordinateXYZM& intPtXYZM = li.getIntersection(i);
+ CoordinateXY intPt(intPtXYZM.x, intPtXYZM.y);
+ /**
+ * Ensure endpoint intersections are added once only, for their canonical segments.
+ * Proper intersections lie on a unique segment so do not need to be checked.
+ * And it is important that the Containing Segment check not be used,
+ * since due to intersection computation roundoff,
+ * it is not reliable in that situation.
+ */
+ if (li.isProper() ||
+ (ssA->isContainingSegment(segIndexA, &intPt) &&
+ ssB->isContainingSegment(segIndexB, &intPt)))
+ {
+ NodeSection* nsa = ssA->createNodeSection(segIndexA, intPt);
+ NodeSection* nsb = ssB->createNodeSection(segIndexB, intPt);
+ topoComputer.addIntersection(nsa, nsb);
+ }
+ }
+}
+
+
+
+} // namespace geos.operation.overlayng
+} // namespace geos.operation
+} // namespace geos
+
+
diff --git a/src/operation/relateng/EdgeSegmentOverlapAction.cpp b/src/operation/relateng/EdgeSegmentOverlapAction.cpp
new file mode 100644
index 000000000..cd32bca78
--- /dev/null
+++ b/src/operation/relateng/EdgeSegmentOverlapAction.cpp
@@ -0,0 +1,53 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2001-2002 Vivid Solutions Inc.
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************
+ *
+ * Last port: index/chain/MonotoneChainOverlapAction.java rev. 1.6 (JTS-1.10)
+ *
+ **********************************************************************/
+
+#include <geos/index/chain/MonotoneChainOverlapAction.h>
+#include <geos/index/chain/MonotoneChain.h>
+#include <geos/noding/SegmentString.h>
+#include <geos/noding/SegmentIntersector.h>
+#include <geos/operation/relateng/EdgeSegmentOverlapAction.h>
+#include <geos/geom/Envelope.h>
+#include <geos/geom/LineSegment.h>
+#include <geos/profiler.h>
+
+//#include <stdio.h>
+
+using geos::index::chain::MonotoneChain;
+using geos::noding::SegmentString;
+using geos::noding::SegmentIntersector;
+
+namespace geos {
+namespace operation { // geos.operation
+namespace relateng { // geos.operation.relateng
+
+
+/* public override */
+void
+EdgeSegmentOverlapAction::overlap(
+ const MonotoneChain& mc1, std::size_t start1,
+ const MonotoneChain& mc2, std::size_t start2)
+{
+ SegmentString* ss1 = static_cast<SegmentString*>(mc1.getContext());
+ SegmentString* ss2 = static_cast<SegmentString*>(mc2.getContext());
+ si.processIntersections(ss1, start1, ss2, start2);
+}
+
+
+} // namespace geos.index.chain
+} // namespace geos.index
+} // namespace geos
diff --git a/src/operation/relateng/EdgeSetIntersector.cpp b/src/operation/relateng/EdgeSetIntersector.cpp
new file mode 100644
index 000000000..826ee94ca
--- /dev/null
+++ b/src/operation/relateng/EdgeSetIntersector.cpp
@@ -0,0 +1,99 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (c) 2024 Martin Davis
+ * Copyright (C) 2024 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#include <geos/geom/Envelope.h>
+#include <geos/geom/Geometry.h>
+#include <geos/noding/SegmentString.h>
+#include <geos/operation/relateng/EdgeSegmentIntersector.h>
+#include <geos/operation/relateng/EdgeSetIntersector.h>
+#include <geos/operation/relateng/EdgeSegmentOverlapAction.h>
+#include <geos/operation/relateng/RelateSegmentString.h>
+#include <geos/index/chain/MonotoneChain.h>
+#include <geos/index/chain/MonotoneChainBuilder.h>
+
+
+using geos::geom::Geometry;
+using geos::geom::Envelope;
+using geos::noding::SegmentString;
+using geos::index::chain::MonotoneChain;
+using geos::index::chain::MonotoneChainBuilder;
+
+
+namespace geos { // geos
+namespace operation { // geos.operation
+namespace relateng { // geos.operation.relateng
+
+
+/* private */
+void
+EdgeSetIntersector::addEdges(std::vector<std::unique_ptr<RelateSegmentString>>& segStrings)
+{
+ for (auto& ss : segStrings) {
+ addToIndex(ss);
+ }
+}
+
+/* private */
+void
+EdgeSetIntersector::addToIndex(std::unique_ptr<RelateSegmentString>& segStr)
+{
+ std::vector<MonotoneChain> segChains;
+ MonotoneChainBuilder::getChains(segStr->getCoordinates(), segStr.get(), segChains);
+
+ for (MonotoneChain& mc : segChains) {
+ if (envelope == nullptr || envelope->intersects(mc.getEnvelope())) {
+ mc.setId(idCounter++);
+ monoChains.push_back(mc);
+ MonotoneChain* mcPtr = &(monoChains.back());
+ index.insert(mcPtr->getEnvelope(), mcPtr);
+ }
+ }
+}
+
+/* public */
+void
+EdgeSetIntersector::process(EdgeSegmentIntersector& intersector)
+{
+ EdgeSegmentOverlapAction overlapAction(intersector);
+
+ for (const MonotoneChain& queryChain : monoChains) {
+
+ std::vector<const MonotoneChain*> overlapChains;
+ index.query(queryChain.getEnvelope(), [&overlapChains](const MonotoneChain* mc) {
+ overlapChains.push_back(mc);
+ });
+
+ for (const MonotoneChain* testChain : overlapChains) {
+ /**
+ * following test makes sure we only compare each pair of chains once
+ * and that we don't compare a chain to itself
+ */
+ if (testChain->getId() <= queryChain.getId())
+ continue;
+
+ testChain->computeOverlaps(&queryChain, &overlapAction);
+ if (intersector.isDone())
+ return;
+ }
+ }
+}
+
+
+
+} // namespace geos.operation.overlayng
+} // namespace geos.operation
+} // namespace geos
+
+
diff --git a/src/operation/relateng/NodeSection.cpp b/src/operation/relateng/NodeSection.cpp
index 81d7d304e..d8d48f544 100644
--- a/src/operation/relateng/NodeSection.cpp
+++ b/src/operation/relateng/NodeSection.cpp
@@ -42,7 +42,7 @@ NodeSection::getVertex(int i) const
/* public */
-const CoordinateXY *
+const CoordinateXY &
NodeSection::nodePt() const
{
return m_nodePt;
@@ -97,6 +97,14 @@ NodeSection::isArea() const
}
+/* public static */
+bool
+NodeSection::isAreaArea(const NodeSection& a, const NodeSection& b)
+{
+ return a.dimension() == Dimension::A && b.dimension() == Dimension::A;
+}
+
+
/* public */
bool
NodeSection::isA() const
@@ -223,9 +231,9 @@ NodeSection::toString() const
if (m_id >= 0) {
ss << "[" << m_id << ":" << m_ringId << "]";
}
- ss << ": " << edgeRep(m_v0, m_nodePt);
+ ss << ": " << edgeRep(m_v0, &m_nodePt);
ss << (m_isNodeAtVertex ? "-V-" : "---");
- ss << " " << edgeRep(m_nodePt, m_v1);
+ ss << " " << edgeRep(&m_nodePt, m_v1);
return ss.str();
}
diff --git a/src/operation/relateng/PolygonNodeConverter.cpp b/src/operation/relateng/PolygonNodeConverter.cpp
index e9b1e9c41..edc582c32 100644
--- a/src/operation/relateng/PolygonNodeConverter.cpp
+++ b/src/operation/relateng/PolygonNodeConverter.cpp
@@ -44,7 +44,7 @@ PolygonNodeConverter::convert(std::vector<const NodeSection*>& polySections)
const NodeSection* ns2)
{
int comp = PolygonNodeTopology::compareAngle(
- ns1->nodePt(),
+ &(ns1->nodePt()),
ns1->getVertex(0),
ns2->getVertex(0));
return comp < 0;
diff --git a/src/operation/relateng/RelateSegmentString.cpp b/src/operation/relateng/RelateSegmentString.cpp
index ccdf46d02..d6f0815a9 100644
--- a/src/operation/relateng/RelateSegmentString.cpp
+++ b/src/operation/relateng/RelateSegmentString.cpp
@@ -99,14 +99,14 @@ RelateSegmentString::getPolygonal() const
/* public */
NodeSection*
-RelateSegmentString::createNodeSection(std::size_t segIndex, const CoordinateXY& intPt) const
+RelateSegmentString::createNodeSection(std::size_t segIndex, const CoordinateXY intPt) const
{
const CoordinateXY& c0 = getCoordinate(segIndex);
const CoordinateXY& c1 = getCoordinate(segIndex + 1);
bool isNodeAtVertex = intPt.equals2D(c0) || intPt.equals2D(c1);
const CoordinateXY* prev = prevVertex(segIndex, &intPt);
const CoordinateXY* next = nextVertex(segIndex, &intPt);
- NodeSection* a = new NodeSection(m_isA, m_dimension, m_id, m_ringId, m_parentPolygonal, isNodeAtVertex, prev, &intPt, next);
+ NodeSection* a = new NodeSection(m_isA, m_dimension, m_id, m_ringId, m_parentPolygonal, isNodeAtVertex, prev, intPt, next);
return a;
}
diff --git a/src/operation/relateng/TopologyComputer.cpp b/src/operation/relateng/TopologyComputer.cpp
new file mode 100644
index 000000000..37adebd6c
--- /dev/null
+++ b/src/operation/relateng/TopologyComputer.cpp
@@ -0,0 +1,546 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (c) 2024 Martin Davis
+ * Copyright (C) 2024 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#include <geos/algorithm/PolygonNodeTopology.h>
+#include <geos/geom/Dimension.h>
+#include <geos/geom/Location.h>
+#include <geos/geom/Position.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/Geometry.h>
+#include <geos/operation/relateng/TopologyComputer.h>
+#include <geos/operation/relateng/TopologyPredicate.h>
+#include <geos/operation/relateng/RelateGeometry.h>
+#include <geos/operation/relateng/RelateNode.h>
+#include <geos/operation/relateng/NodeSection.h>
+#include <geos/util/IllegalStateException.h>
+#include <sstream>
+
+
+using geos::algorithm::PolygonNodeTopology;
+using geos::geom::Coordinate;
+using geos::geom::CoordinateXY;
+using geos::geom::Geometry;
+using geos::geom::Dimension;
+using geos::geom::Location;
+using geos::geom::Position;
+using geos::util::IllegalStateException;
+
+
+namespace geos { // geos
+namespace operation { // geos.operation
+namespace relateng { // geos.operation.relateng
+
+
+/* private */
+void
+TopologyComputer::initExteriorDims()
+{
+ int dimRealA = geomA->getDimensionReal();
+ int dimRealB = geomB->getDimensionReal();
+
+ /**
+ * For P/L case, P exterior intersects L interior
+ */
+ if (dimRealA == Dimension::P && dimRealB == Dimension::L) {
+ updateDim(Location::EXTERIOR, Location::INTERIOR, Dimension::L);
+ }
+ else if (dimRealA == Dimension::L && dimRealB == Dimension::P) {
+ updateDim(Location::INTERIOR, Location::EXTERIOR, Dimension::L);
+ }
+ /**
+ * For P/A case, the Area Int and Bdy intersect the Point exterior.
+ */
+ else if (dimRealA == Dimension::P && dimRealB == Dimension::A) {
+ updateDim(Location::EXTERIOR, Location::INTERIOR, Dimension::A);
+ updateDim(Location::EXTERIOR, Location::BOUNDARY, Dimension::L);
+ }
+ else if (dimRealA == Dimension::A && dimRealB == Dimension::P) {
+ updateDim(Location::INTERIOR, Location::EXTERIOR, Dimension::A);
+ updateDim(Location::BOUNDARY, Location::EXTERIOR, Dimension::L);
+ }
+ else if (dimRealA == Dimension::L && dimRealB == Dimension::A) {
+ updateDim(Location::EXTERIOR, Location::INTERIOR, Dimension::A);
+ }
+ else if (dimRealA == Dimension::A && dimRealB == Dimension::L) {
+ updateDim(Location::INTERIOR, Location::EXTERIOR, Dimension::A);
+ }
+ //-- cases where one geom is EMPTY
+ else if (dimRealA == Dimension::False || dimRealB == Dimension::False) {
+ if (dimRealA != Dimension::False) {
+ initExteriorEmpty(RelateGeometry::GEOM_A);
+ }
+ if (dimRealB != Dimension::False) {
+ initExteriorEmpty(RelateGeometry::GEOM_B);
+ }
+ }
+}
+
+
+/* private */
+void
+TopologyComputer::initExteriorEmpty(bool geomNonEmpty)
+{
+ int dimNonEmpty = getDimension(geomNonEmpty);
+ switch (dimNonEmpty) {
+ case Dimension::P:
+ updateDim(geomNonEmpty, Location::INTERIOR, Location::EXTERIOR, Dimension::P);
+ break;
+ case Dimension::L:
+ if (getGeometry(geomNonEmpty)->hasBoundary()) {
+ updateDim(geomNonEmpty, Location::BOUNDARY, Location::EXTERIOR, Dimension::P);
+ }
+ updateDim(geomNonEmpty, Location::INTERIOR, Location::EXTERIOR, Dimension::L);
+ break;
+ case Dimension::A:
+ updateDim(geomNonEmpty, Location::BOUNDARY, Location::EXTERIOR, Dimension::L);
+ updateDim(geomNonEmpty, Location::INTERIOR, Location::EXTERIOR, Dimension::A);
+ break;
+ }
+}
+
+
+/* private */
+RelateGeometry*
+TopologyComputer::getGeometry(bool isA) const
+{
+ return isA ? geomA : geomB;
+}
+
+
+/* public */
+int
+TopologyComputer::getDimension(bool isA) const
+{
+ return getGeometry(isA)->getDimension();
+}
+
+
+/* public */
+bool
+TopologyComputer::isAreaArea() const
+{
+ return getDimension(RelateGeometry::GEOM_A) == Dimension::A
+ && getDimension(RelateGeometry::GEOM_B) == Dimension::A;
+}
+
+
+/* public */
+bool
+TopologyComputer::isSelfNodingRequired() const
+{
+ //TODO: change to testing for lines or GC with > 1 polygon
+ if (geomA->isPointsOrPolygons()) return false;
+ if (geomB->isPointsOrPolygons()) return false;
+ return predicate->requireSelfNoding();
+}
+
+
+/* public */
+bool
+TopologyComputer::isExteriorCheckRequired(bool isA) const
+{
+ return predicate->requireExteriorCheck(isA);
+}
+
+
+/* private */
+void
+TopologyComputer::updateDim(Location locA, Location locB, int dimension)
+{
+ predicate->updateDimension(locA, locB, dimension);
+}
+
+
+/* private */
+void
+TopologyComputer::updateDim(bool isAB, Location loc1, Location loc2, int dimension)
+{
+ if (isAB) {
+ updateDim(loc1, loc2, dimension);
+ }
+ else {
+ // is ordered BA
+ updateDim(loc2, loc1, dimension);
+ }
+}
+
+
+/* public */
+bool
+TopologyComputer::isResultKnown() const
+{
+ return predicate->isKnown();
+}
+
+
+/* public */
+bool
+TopologyComputer::getResult() const
+{
+ return predicate->value();
+}
+
+
+/* public */
+void
+TopologyComputer::finish()
+{
+ predicate->finish();
+}
+
+
+/* private */
+NodeSections *
+TopologyComputer::getNodeSections(const CoordinateXY& nodePt)
+{
+ NodeSections* ns;
+ auto result = nodeMap.find(nodePt);
+ if (result == nodeMap.end()) {
+ ns = new NodeSections(&nodePt);
+ nodeSectionsStore.emplace_back(ns);
+ nodeMap[nodePt] = ns;
+ }
+ else {
+ ns = result->second;
+ }
+ return ns;
+}
+
+
+/* public */
+void
+TopologyComputer::addIntersection(NodeSection* a, NodeSection* b)
+{
+ if (! a->isSameGeometry(b)) {
+ updateIntersectionAB(a, b);
+ }
+ //-- add edges to node to allow full topology evaluation later
+ addNodeSections(a, b);
+}
+
+
+/* private */
+void
+TopologyComputer::updateIntersectionAB(const NodeSection* a, const NodeSection* b)
+{
+ if (NodeSection::isAreaArea(*a, *b)) {
+ updateAreaAreaCross(a, b);
+ }
+ updateNodeLocation(a, b);
+}
+
+
+/* private */
+void
+TopologyComputer::updateAreaAreaCross(const NodeSection* a, const NodeSection* b)
+{
+ bool isProper = NodeSection::isProper(*a, *b);
+ if (isProper || PolygonNodeTopology::isCrossing(&(a->nodePt()),
+ a->getVertex(0), a->getVertex(1),
+ b->getVertex(0), b->getVertex(1)))
+ {
+ updateDim(Location::INTERIOR, Location::INTERIOR, Dimension::A);
+ }
+}
+
+
+/* private */
+void
+TopologyComputer::updateNodeLocation(const NodeSection* a, const NodeSection* b)
+{
+ const CoordinateXY& pt = a->nodePt();
+ Location locA = geomA->locateNode(&pt, a->getPolygonal());
+ Location locB = geomB->locateNode(&pt, b->getPolygonal());
+ updateDim(locA, locB, Dimension::P);
+}
+
+/* private */
+void
+TopologyComputer::addNodeSections(NodeSection* ns0, NodeSection* ns1)
+{
+ NodeSections* sections = getNodeSections(ns0->nodePt());
+ sections->addNodeSection(ns0);
+ sections->addNodeSection(ns1);
+}
+
+/* public */
+void
+TopologyComputer::addPointOnPointInterior()
+{
+ updateDim(Location::INTERIOR, Location::INTERIOR, Dimension::P);
+}
+
+/* public */
+void
+TopologyComputer::addPointOnPointExterior(bool isGeomA)
+{
+ updateDim(isGeomA, Location::INTERIOR, Location::EXTERIOR, Dimension::P);
+}
+
+/* public */
+void
+TopologyComputer::addPointOnGeometry(bool isA, Location locTarget, int dimTarget)
+{
+ updateDim(isA, Location::INTERIOR, locTarget, Dimension::P);
+ switch (dimTarget) {
+ case Dimension::P:
+ return;
+ case Dimension::L:
+ /**
+ * Because zero-length lines are handled,
+ * a point lying in the exterior of the line target
+ * may imply either P or L for the Exterior interaction
+ */
+ //TODO: determine if effective dimension of linear target is L?
+ //updateDim(isGeomA, Location::EXTERIOR, locTarget, Dimension::P);
+ return;
+ case Dimension::A:
+ /**
+ * If a point intersects an area target, then the area interior and boundary
+ * must extend beyond the point and thus interact with its exterior.
+ */
+ updateDim(isA, Location::EXTERIOR, Location::INTERIOR, Dimension::A);
+ updateDim(isA, Location::EXTERIOR, Location::BOUNDARY, Dimension::L);
+ return;
+ }
+ throw IllegalStateException("Unknown target dimension: " + std::to_string(dimTarget));
+}
+
+/* public */
+void
+TopologyComputer::addLineEndOnGeometry(bool isLineA, Location locLineEnd, Location locTarget, int dimTarget)
+{
+ switch (dimTarget) {
+ case Dimension::P:
+ addLineEndOnPoint(isLineA, locLineEnd, locTarget);
+ return;
+ case Dimension::L:
+ addLineEndOnLine(isLineA, locLineEnd, locTarget);
+ return;
+ case Dimension::A:
+ addLineEndOnArea(isLineA, locLineEnd, locTarget);
+ return;
+ }
+ throw IllegalStateException("Unknown target dimension: " + std::to_string(dimTarget));
+}
+
+/* private */
+void
+TopologyComputer::addLineEndOnPoint(bool isLineA, Location locLineEnd, Location locPoint)
+{
+ updateDim(isLineA, locLineEnd, locPoint, Dimension::P);
+}
+
+/* private */
+void
+TopologyComputer::addLineEndOnLine(bool isLineA, Location locLineEnd, Location locLine)
+{
+ updateDim(isLineA, locLineEnd, locLine, Dimension::P);
+ /**
+ * When a line end is in the exterior, some length of the line interior
+ * must also be in the exterior.
+ * This works for zero-length lines as well.
+ */
+
+ if (locLine == Location::EXTERIOR) {
+ updateDim(isLineA, Location::INTERIOR, Location::EXTERIOR, Dimension::L);
+ }
+}
+
+/* private */
+void
+TopologyComputer::addLineEndOnArea(bool isLineA, Location locLineEnd, Location locArea)
+{
+ if (locArea == Location::BOUNDARY) {
+ updateDim(isLineA, locLineEnd, locArea, Dimension::P);
+ }
+ else {
+ //TODO: handle zero-length lines?
+ updateDim(isLineA, Location::INTERIOR, locArea, Dimension::L);
+ updateDim(isLineA, locLineEnd, locArea, Dimension::P);
+ updateDim(isLineA, Location::EXTERIOR, locArea, Dimension::A);
+ }
+}
+
+
+/* public */
+void
+TopologyComputer::addAreaVertex(bool isAreaA, Location locArea, Location locTarget, int dimTarget)
+{
+ if (locTarget == Location::EXTERIOR) {
+ updateDim(isAreaA, Location::INTERIOR, Location::EXTERIOR, Dimension::A);
+ /**
+ * If area vertex is on Boundary further topology can be deduced
+ * from the neighbourhood around the boundary vertex.
+ * This is always the case for polygonal geometries.
+ * For GCs, the vertex may be either on boundary or in interior
+ * (i.e. of overlapping or adjacent polygons)
+ */
+ if (locArea == Location::BOUNDARY) {
+ updateDim(isAreaA, Location::BOUNDARY, Location::EXTERIOR, Dimension::L);
+ updateDim(isAreaA, Location::EXTERIOR, Location::EXTERIOR, Dimension::A);
+ }
+ return;
+ }
+
+ switch (dimTarget) {
+ case Dimension::P:
+ addAreaVertexOnPoint(isAreaA, locArea);
+ return;
+ case Dimension::L:
+ addAreaVertexOnLine(isAreaA, locArea, locTarget);
+ return;
+ case Dimension::A:
+ addAreaVertexOnArea(isAreaA, locArea, locTarget);
+ return;
+ }
+ throw IllegalStateException("Unknown target dimension: " + std::to_string(dimTarget));
+}
+
+
+/* private */
+void
+TopologyComputer::addAreaVertexOnPoint(bool isAreaA, Location locArea)
+{
+ //-- Assert: locArea != EXTERIOR
+ //-- Assert: locTarget == INTERIOR
+ /**
+ * The vertex location intersects the Point.
+ */
+ updateDim(isAreaA, locArea, Location::INTERIOR, Dimension::P);
+ /**
+ * The area interior intersects the point's exterior neighbourhood.
+ */
+ updateDim(isAreaA, Location::INTERIOR, Location::EXTERIOR, Dimension::A);
+ /**
+ * If the area vertex is on the boundary,
+ * the area boundary and exterior intersect the point's exterior neighbourhood
+ */
+ if (locArea == Location::BOUNDARY) {
+ updateDim(isAreaA, Location::BOUNDARY, Location::EXTERIOR, Dimension::L);
+ updateDim(isAreaA, Location::EXTERIOR, Location::EXTERIOR, Dimension::A);
+ }
+}
+
+
+/* private */
+void
+TopologyComputer::addAreaVertexOnLine(bool isAreaA, Location locArea, Location locTarget)
+{
+ //-- Assert: locArea != EXTERIOR
+ /**
+ * If an area vertex intersects a line, all we know is the
+ * intersection at that point.
+ * e.g. the line may or may not be collinear with the area boundary,
+ * and the line may or may not intersect the area interior.
+ * Full topology is determined later by node analysis
+ */
+ updateDim(isAreaA, locArea, locTarget, Dimension::P);
+ if (locArea == Location::INTERIOR) {
+ /**
+ * The area interior intersects the line's exterior neighbourhood.
+ */
+ updateDim(isAreaA, Location::INTERIOR, Location::EXTERIOR, Dimension::A);
+ }
+}
+
+
+/* public */
+void
+TopologyComputer::addAreaVertexOnArea(bool isAreaA, Location locArea, Location locTarget)
+{
+ if (locTarget == Location::BOUNDARY) {
+ if (locArea == Location::BOUNDARY) {
+ //-- B/B topology is fully computed later by node analysis
+ updateDim(isAreaA, Location::BOUNDARY, Location::BOUNDARY, Dimension::P);
+ }
+ else {
+ // locArea == INTERIOR
+ updateDim(isAreaA, Location::INTERIOR, Location::INTERIOR, Dimension::A);
+ updateDim(isAreaA, Location::INTERIOR, Location::BOUNDARY, Dimension::L);
+ updateDim(isAreaA, Location::INTERIOR, Location::EXTERIOR, Dimension::A);
+ }
+ }
+ else {
+ //-- locTarget is INTERIOR or EXTERIOR`
+ updateDim(isAreaA, Location::INTERIOR, locTarget, Dimension::A);
+ /**
+ * If area vertex is on Boundary further topology can be deduced
+ * from the neighbourhood around the boundary vertex.
+ * This is always the case for polygonal geometries.
+ * For GCs, the vertex may be either on boundary or in interior
+ * (i.e. of overlapping or adjacent polygons)
+ */
+ if (locArea == Location::BOUNDARY) {
+ updateDim(isAreaA, Location::BOUNDARY, locTarget, Dimension::L);
+ updateDim(isAreaA, Location::EXTERIOR, locTarget, Dimension::A);
+ }
+ }
+}
+
+
+/* public */
+void
+TopologyComputer::evaluateNodes()
+{
+ for (auto& kv : nodeMap) {
+ NodeSections* nodeSections = kv.second;
+ if (nodeSections->hasInteractionAB()) {
+ evaluateNode(nodeSections);
+ if (isResultKnown())
+ return;
+ }
+ }
+}
+
+
+/* private */
+void
+TopologyComputer::evaluateNode(NodeSections* nodeSections)
+{
+ const CoordinateXY* p = nodeSections->getCoordinate();
+ std::unique_ptr<RelateNode> node = nodeSections->createNode();
+ //-- Node must have edges for geom, but may also be in interior of a overlapping GC
+ bool isAreaInteriorA = geomA->isNodeInArea(p, nodeSections->getPolygonal(RelateGeometry::GEOM_A));
+ bool isAreaInteriorB = geomB->isNodeInArea(p, nodeSections->getPolygonal(RelateGeometry::GEOM_B));
+ node->finish(isAreaInteriorA, isAreaInteriorB);
+ evaluateNodeEdges(node.get());
+}
+
+
+/* private */
+void
+TopologyComputer::evaluateNodeEdges(const RelateNode* node)
+{
+ //TODO: collect distinct dim settings by using temporary matrix?
+ for (const std::unique_ptr<RelateEdge>& e : node->getEdges()) {
+ //-- An optimization to avoid updates for cases with a linear geometry
+ if (isAreaArea()) {
+ updateDim(e->location(RelateGeometry::GEOM_A, Position::LEFT),
+ e->location(RelateGeometry::GEOM_B, Position::LEFT), Dimension::A);
+ updateDim(e->location(RelateGeometry::GEOM_A, Position::RIGHT),
+ e->location(RelateGeometry::GEOM_B, Position::RIGHT), Dimension::A);
+ }
+ updateDim(e->location(RelateGeometry::GEOM_A, Position::ON),
+ e->location(RelateGeometry::GEOM_B, Position::ON), Dimension::L);
+ }
+}
+
+
+} // namespace geos.operation.overlayng
+} // namespace geos.operation
+} // namespace geos
+
+
diff --git a/tests/unit/operation/relateng/PolygonNodeConverterTest.cpp b/tests/unit/operation/relateng/PolygonNodeConverterTest.cpp
index 6037df049..2d3ee22d7 100644
--- a/tests/unit/operation/relateng/PolygonNodeConverterTest.cpp
+++ b/tests/unit/operation/relateng/PolygonNodeConverterTest.cpp
@@ -41,7 +41,7 @@ struct test_polygonnodeconverter_data {
NodeSection* section(int ringId, double v0x, double v0y, double nx, double ny, double v1x, double v1y) {
return new NodeSection(true, Dimension::A, 1, ringId, nullptr, false,
- new Coordinate(v0x, v0y), new Coordinate(nx, ny), new Coordinate(v1x, v1y));
+ new Coordinate(v0x, v0y), Coordinate(nx, ny), new Coordinate(v1x, v1y));
}
std::vector<const NodeSection*>
@@ -85,7 +85,6 @@ struct test_polygonnodeconverter_data {
for (std::size_t i = 0; i < ns.size(); i++) {
delete ns[i]->getVertex(0);
delete ns[i]->getVertex(1);
- delete ns[i]->nodePt();
delete ns[i];
}
}
diff --git a/tests/unit/operation/relateng/RelateGeometryTest.cpp b/tests/unit/operation/relateng/RelateGeometryTest.cpp
index b3718ec3a..8ecd91e1f 100644
--- a/tests/unit/operation/relateng/RelateGeometryTest.cpp
+++ b/tests/unit/operation/relateng/RelateGeometryTest.cpp
@@ -123,7 +123,12 @@ void object::test<9> ()
checkDimension("GEOMETRYCOLLECTION (POLYGON EMPTY, LINESTRING (1 1, 5 4), POINT (6 5))", 2, 1);
}
-
+template<>
+template<>
+void object::test<10> ()
+{
+ checkDimension("LINESTRING (0 0, 0 0, 9 9)", 1, 1);
+}
-----------------------------------------------------------------------
Summary of changes:
include/geos/index/chain/MonotoneChain.h | 7 +
.../operation/relateng/EdgeSegmentIntersector.h | 80 +++
.../operation/relateng/EdgeSegmentOverlapAction.h | 75 +++
.../geos/operation/relateng/EdgeSetIntersector.h | 95 ++++
include/geos/operation/relateng/NodeSection.h | 13 +-
.../geos/operation/relateng/RelateSegmentString.h | 2 +-
include/geos/operation/relateng/TopologyComputer.h | 226 +++++++++
src/operation/relateng/AdjacentEdgeLocator.cpp | 2 +-
src/operation/relateng/DimensionLocation.cpp | 8 +-
src/operation/relateng/EdgeSegmentIntersector.cpp | 114 +++++
.../relateng/EdgeSegmentOverlapAction.cpp} | 25 +-
src/operation/relateng/EdgeSetIntersector.cpp | 99 ++++
src/operation/relateng/NodeSection.cpp | 14 +-
src/operation/relateng/PolygonNodeConverter.cpp | 2 +-
src/operation/relateng/RelateSegmentString.cpp | 4 +-
src/operation/relateng/TopologyComputer.cpp | 546 +++++++++++++++++++++
.../relateng/PolygonNodeConverterTest.cpp | 3 +-
.../unit/operation/relateng/RelateGeometryTest.cpp | 7 +-
18 files changed, 1293 insertions(+), 29 deletions(-)
create mode 100644 include/geos/operation/relateng/EdgeSegmentIntersector.h
create mode 100644 include/geos/operation/relateng/EdgeSegmentOverlapAction.h
create mode 100644 include/geos/operation/relateng/EdgeSetIntersector.h
create mode 100644 include/geos/operation/relateng/TopologyComputer.h
create mode 100644 src/operation/relateng/EdgeSegmentIntersector.cpp
copy src/{index/chain/MonotoneChainOverlapAction.cpp => operation/relateng/EdgeSegmentOverlapAction.cpp} (57%)
create mode 100644 src/operation/relateng/EdgeSetIntersector.cpp
create mode 100644 src/operation/relateng/TopologyComputer.cpp
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list