[geos-commits] [SCM] GEOS branch main-relate-ng updated. 4f213c1df393628e531622b49337e03184869e27

git at osgeo.org git at osgeo.org
Sun Jul 28 16:45:19 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  4f213c1df393628e531622b49337e03184869e27 (commit)
       via  5b8e2679e2b0db5d458007abe392ee0a8d8fd8d2 (commit)
       via  4907d593aa6aa31c25dbe73933bf66995a04e498 (commit)
      from  8ee1ae84dfdfe26101b74271fd3c426c8593df56 (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 4f213c1df393628e531622b49337e03184869e27
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date:   Sun Jul 28 16:44:38 2024 -0700

    Add RelateNG and RelateNGTest

diff --git a/include/geos/geom/GeometryCollection.h b/include/geos/geom/GeometryCollection.h
index 27de4c7cf..9437f6af3 100644
--- a/include/geos/geom/GeometryCollection.h
+++ b/include/geos/geom/GeometryCollection.h
@@ -196,13 +196,6 @@ public:
         return &envelope;
     }
 
-    /**
-     * \brief
-     * Recurse into collection and populate vector with just the
-     * simple non-collection components of the collection.
-     */
-    void getAllGeometries(std::vector<const Geometry*>& geoms) const;
-
 protected:
 
     GeometryCollection(const GeometryCollection& gc);
diff --git a/include/geos/geom/util/GeometryLister.h b/include/geos/geom/util/GeometryLister.h
new file mode 100644
index 000000000..912518ef1
--- /dev/null
+++ b/include/geos/geom/util/GeometryLister.h
@@ -0,0 +1,92 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2011 Sandro Santilli <strk at kbt.io>
+ * Copyright (C) 2001-2002 Vivid Solutions Inc.
+ * Copyright (C) 2006 Refractions Research 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: geom/util/GeometryExtracter.java r320 (JTS-1.12)
+ *
+ **********************************************************************/
+
+#pragma once
+
+#include <geos/export.h>
+#include <geos/geom/GeometryFilter.h>
+#include <geos/geom/GeometryCollection.h>
+#include <vector>
+
+namespace geos {
+namespace geom { // geos.geom
+namespace util { // geos.geom.util
+
+/**
+ * Extracts all the components of a collection, or just echoes back a
+ * pointers to singletons.
+ */
+class GEOS_DLL GeometryLister {
+
+public:
+
+    /**
+     * Extracts the components from a {@link Geometry}
+     * and adds them to the provided container.
+     *
+     * Useful for iterating over components of a collection.
+     *
+     * @param geom the geometry from which to extract
+     * @param lst the list to add the extracted elements to
+     */
+    static void
+    list(const Geometry* geom, std::vector<const Geometry*>& lst)
+    {
+        if(geom->isCollection()) {
+            GeometryLister::Lister lister(lst);
+            geom->apply_ro(&lister);
+        }
+        else {
+            lst.push_back(geom);
+        }
+    }
+
+private:
+
+    struct Lister : public GeometryFilter {
+
+        /**
+         * Constructs a filter with a list in which to store the elements found.
+         *
+         * @param comps the container to extract into (will push_back to it)
+         */
+        Lister(std::vector<const Geometry*>& p_geoms) : geoms(p_geoms) {}
+
+        std::vector<const Geometry*>& geoms;
+
+        void
+        filter_ro(const Geometry* geom) override
+        {
+            if(!geom->isCollection()) {
+                geoms.push_back(geom);
+            }
+        }
+
+        // // Declare type as noncopyable
+        // Lister(const Lister& other);
+        // Lister& operator=(const Lister& rhs);
+    };
+};
+
+
+} // namespace geos.geom.util
+} // namespace geos.geom
+} // namespace geos
+
diff --git a/include/geos/operation/relateng/EdgeSegmentIntersector.h b/include/geos/operation/relateng/EdgeSegmentIntersector.h
index 5153b68fd..1b05b6f13 100644
--- a/include/geos/operation/relateng/EdgeSegmentIntersector.h
+++ b/include/geos/operation/relateng/EdgeSegmentIntersector.h
@@ -62,8 +62,8 @@ private:
 
 public:
 
-    EdgeSegmentIntersector(TopologyComputer& p_topoBuilder)
-        : topoComputer(p_topoBuilder)
+    EdgeSegmentIntersector(TopologyComputer& p_topoComputer)
+        : topoComputer(p_topoComputer)
         {};
 
     void processIntersections(
diff --git a/include/geos/operation/relateng/EdgeSetIntersector.h b/include/geos/operation/relateng/EdgeSetIntersector.h
index bf95bb6d7..70b6cb315 100644
--- a/include/geos/operation/relateng/EdgeSetIntersector.h
+++ b/include/geos/operation/relateng/EdgeSetIntersector.h
@@ -64,16 +64,16 @@ private:
 
     // Methods
 
-    void addToIndex(std::unique_ptr<RelateSegmentString>& segStr);
+    void addToIndex(const SegmentString* segStr);
 
-    void addEdges(std::vector<std::unique_ptr<RelateSegmentString>>& segStrings);
+    void addEdges(std::vector<const SegmentString*>& segStrings);
 
 
 public:
 
     EdgeSetIntersector(
-        std::vector<std::unique_ptr<RelateSegmentString>>& edgesA,
-        std::vector<std::unique_ptr<RelateSegmentString>>& edgesB,
+        std::vector<const SegmentString*>& edgesA,
+        std::vector<const SegmentString*>& edgesB,
         const Envelope* env)
         : envelope(env)
         , idCounter(0)
diff --git a/include/geos/operation/relateng/IntersectionMatrixPattern.h b/include/geos/operation/relateng/IntersectionMatrixPattern.h
index a7ef3462c..9e7f2f504 100644
--- a/include/geos/operation/relateng/IntersectionMatrixPattern.h
+++ b/include/geos/operation/relateng/IntersectionMatrixPattern.h
@@ -41,20 +41,20 @@ public:
      * A DE-9IM pattern to detect whether two polygonal geometries are adjacent along
      * an edge, but do not overlap.
      */
-    static constexpr std::string ADJACENT = "F***1****";
+    static constexpr const char* ADJACENT = "F***1****";
 
     /**
      * A DE-9IM pattern to detect a geometry which properly contains another
      * geometry (i.e. which lies entirely in the interior of the first geometry).
      */
-    static constexpr std::string CONTAINS_PROPERLY = "T**FF*FF*";
+    static constexpr const char* CONTAINS_PROPERLY = "T**FF*FF*";
 
     /**
      * A DE-9IM pattern to detect if two geometries intersect in their interiors.
      * This can be used to determine if a polygonal coverage contains any overlaps
      * (although not whether they are correctly noded).
      */
-    static constexpr std::string INTERIOR_INTERSECTS = "T********";
+    static constexpr const char* INTERIOR_INTERSECTS = "T********";
 
 
 };
diff --git a/include/geos/operation/relateng/NodeSections.h b/include/geos/operation/relateng/NodeSections.h
index 9f321dd27..d7d08b342 100644
--- a/include/geos/operation/relateng/NodeSections.h
+++ b/include/geos/operation/relateng/NodeSections.h
@@ -17,7 +17,7 @@
 
 #include <vector>
 #include <memory>
-// #include <geos/operation/relateng/NodeSection.h>
+#include <geos/operation/relateng/NodeSection.h>
 #include <geos/export.h>
 
 
@@ -26,7 +26,7 @@ namespace geos {
 namespace operation {
 namespace relateng {
 class RelateNode;
-class NodeSection;
+// class NodeSection;
 }
 }
 namespace geom {
diff --git a/include/geos/operation/relateng/RelateGeometry.h b/include/geos/operation/relateng/RelateGeometry.h
index 59ebddcbb..ae4101111 100644
--- a/include/geos/operation/relateng/RelateGeometry.h
+++ b/include/geos/operation/relateng/RelateGeometry.h
@@ -20,9 +20,9 @@
 #include <geos/geom/Dimension.h>
 #include <geos/geom/Location.h>
 #include <geos/operation/relateng/RelatePointLocator.h>
+#include <geos/operation/relateng/RelateSegmentString.h>
 #include <geos/export.h>
 
-
 #include <string>
 #include <sstream>
 
@@ -38,15 +38,14 @@ namespace geom {
     class MultiPolygon;
     class Point;
 }
-namespace operation {
-namespace relateng {
-    class RelateSegmentString;
-}
+namespace noding {
+    class SegmentString;
 }
 }
 
 
 using geos::algorithm::BoundaryNodeRule;
+using geos::noding::SegmentString;
 using namespace geos::geom;
 
 
@@ -76,6 +75,12 @@ private:
     bool hasLines = false;
     bool hasAreas = false;
 
+    /*
+     * Memory contexts for lower level allocations
+     */
+    std::vector<std::unique_ptr<const RelateSegmentString>> segStringStore;
+    std::vector<std::unique_ptr<CoordinateSequence>> csStore;
+
 
     // Methods
 
@@ -99,17 +104,22 @@ private:
     void extractSegmentStringsFromAtomic(bool isA,
         const Geometry* geom, const MultiPolygon* parentPolygonal,
         const Envelope* env,
-        std::vector<std::unique_ptr<RelateSegmentString>>& segStrings);
+        std::vector<const SegmentString*>& segStrings);
 
     void extractRingToSegmentString(bool isA,
         const LinearRing* ring, int ringId, const Envelope* env,
         const Geometry* parentPoly,
-        std::vector<std::unique_ptr<RelateSegmentString>>& segStrings);
+        std::vector<const SegmentString*>& segStrings);
 
     void extractSegmentStrings(bool isA,
         const Envelope* env, const Geometry* geom,
-        std::vector<std::unique_ptr<RelateSegmentString>>& segStrings);
+        std::vector<const SegmentString*>& segStrings);
 
+    const CoordinateSequence* orientAndRemoveRepeated(
+        const CoordinateSequence* cs, bool orientCW);
+
+    const CoordinateSequence* removeRepeated(
+        const CoordinateSequence* cs);
 
 public:
 
@@ -191,20 +201,28 @@ public:
     std::vector<const Point*> getEffectivePoints();
 
     /**
-     * Extract RelateSegmentStrings from the geometry which
+     * Extract RSegmentStrings from the geometry which
      * intersect a given envelope.
      * If the envelope is null all edges are extracted.
      * @param geomA
      *
      * @param env the envelope to extract around (may be null)
-     * @return a list of RelateSegmentStrings
+     * @return a list of SegmentStrings
      */
-    std::vector<std::unique_ptr<RelateSegmentString>> extractSegmentStrings(bool isA, const Envelope* env);
+    std::vector<const SegmentString*> extractSegmentStrings(bool isA, const Envelope* env);
 
     std::string toString() const;
 
     friend std::ostream& operator<<(std::ostream& os, const RelateGeometry& rg);
 
+    /**
+     * Disable copy construction and assignment. Needed to make this
+     * class compile under MSVC. (See https://stackoverflow.com/q/29565299)
+     * Classes with members that are vector<> of unique_ptr<> need this.
+     */
+    RelateGeometry(const RelateGeometry&) = delete;
+    RelateGeometry& operator=(const RelateGeometry&) = delete;
+
 };
 
 } // namespace geos.operation.relateng
diff --git a/include/geos/operation/relateng/RelateNG.h b/include/geos/operation/relateng/RelateNG.h
new file mode 100644
index 000000000..60e637739
--- /dev/null
+++ b/include/geos/operation/relateng/RelateNG.h
@@ -0,0 +1,261 @@
+/**********************************************************************
+ *
+ * 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/export.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/noding/MCIndexSegmentSetMutualIntersector.h>
+#include <geos/operation/relateng/RelateGeometry.h>
+#include <string>
+#include <sstream>
+
+
+// Forward declarations
+namespace geos {
+namespace algorithm {
+    class BoundaryNodeRule;
+
+}
+namespace geom {
+    class Geometry;
+}
+namespace noding {
+}
+namespace operation {
+namespace relateng {
+    class TopologyPredicate;
+    class TopologyComputer;
+    class EdgeSegmentIntersector;
+}
+}
+}
+
+
+using geos::geom::CoordinateXY;
+using geos::geom::CoordinateXY;
+using geos::geom::Geometry;
+using geos::algorithm::BoundaryNodeRule;
+using geos::noding::MCIndexSegmentSetMutualIntersector;
+
+
+namespace geos {      // geos.
+namespace operation { // geos.operation
+namespace relateng { // geos.operation.relateng
+
+
+/**
+ * Computes the value of topological predicates between two geometries based on the
+ * Dimensionally-Extended 9-Intersection Model <https://en.wikipedia.org/wiki/DE-9IM> (DE-9IM).
+ * Standard and custom topological predicates are provided by RelatePredicate.
+ *
+ * The RelateNG algorithm has the following capabilities:
+ *
+ *   * Efficient short-circuited evaluation of topological predicates
+ *     (including matching custom DE-9IM matrix patterns)
+ *   * Optimized repeated evaluation of predicates against a single geometry
+ *     via cached spatial indexes (AKA "prepared mode")
+ *   * Robust computation (only point-local topology is required,
+ *     so invalid geometry topology does not cause failures)
+ *   * GeometryCollection inputs containing mixed types and overlapping polygons
+ *     are supported, using union semantics.
+ *   * Zero-length LineStrings are treated as being topologically identical to Points.
+ *   * Support for BoundaryNodeRule.
+ *
+ * See IntersectionMatrixPattern for a description of DE-9IM patterns.
+ *
+ * If not specified, the standard BoundaryNodeRule::MOD2_BOUNDARY_RULE is used.
+ *
+ * RelateNG operates in 2D only; it ignores any Z ordinates.
+ *
+ * This implementation replaces RelateOp and PreparedGeometry.
+ *
+ * FUTURE WORK
+ *
+ *   * Support for a distance tolerance to provide "approximate" predicate evaluation
+ *
+ * @author Martin Davis
+ *
+ * @see RelateOp
+ * @see PreparedGeometry
+ */
+class GEOS_DLL RelateNG {
+
+private:
+
+    // Members
+    const BoundaryNodeRule& boundaryNodeRule;
+    RelateGeometry geomA;
+    std::unique_ptr<MCIndexSegmentSetMutualIntersector> edgeMutualInt = nullptr;
+
+    // Methods
+
+    RelateNG(const Geometry* inputA, bool isPrepared, const BoundaryNodeRule& bnRule)
+        : boundaryNodeRule(bnRule)
+        , geomA(inputA, isPrepared, bnRule)
+        {}
+
+    RelateNG(const Geometry* inputA, bool isPrepared)
+        : RelateNG(inputA, isPrepared, BoundaryNodeRule::getBoundaryRuleMod2())
+        {}
+
+    bool hasRequiredEnvelopeInteraction(const Geometry* b, TopologyPredicate& predicate);
+
+    bool finishValue(TopologyPredicate& predicate);
+
+    void computePP(RelateGeometry& geomB, TopologyComputer& topoComputer);
+
+    void computeAtPoints(RelateGeometry& geom, bool isA, RelateGeometry& geomTarget, TopologyComputer& topoComputer);
+
+    bool computePoints(RelateGeometry& geom, bool isA, RelateGeometry& geomTarget, TopologyComputer& topoComputer);
+
+    void computePoint(bool isA, const CoordinateXY* pt, RelateGeometry& geomTarget, TopologyComputer& topoComputer);
+
+    bool computeLineEnds(RelateGeometry& geom, bool isA, RelateGeometry& geomTarget, TopologyComputer& topoComputer);
+
+    bool computeLineEnd(RelateGeometry& geom, bool isA, const CoordinateXY* pt, RelateGeometry& geomTarget, TopologyComputer& topoComputer);
+
+    bool computeAreaVertex(RelateGeometry& geom, bool isA, RelateGeometry& geomTarget, TopologyComputer& topoComputer);
+
+    bool computeAreaVertex(RelateGeometry& geom, bool isA, const LinearRing* ring, RelateGeometry& geomTarget, TopologyComputer& topoComputer);
+
+    void computeAtEdges(RelateGeometry& geomB, TopologyComputer& topoComputer);
+
+    void computeEdgesAll(std::vector<const SegmentString*>& edgesB, const Envelope* envInt, EdgeSegmentIntersector& intersector);
+
+    void computeEdgesMutual(std::vector<const SegmentString*>& edgesB, const Envelope* envInt, EdgeSegmentIntersector& intersector);
+
+
+
+public:
+
+    /**
+     * Tests whether the topological relationship between two geometries
+     * satisfies a topological predicate.
+     *
+     * @param a the A input geometry
+     * @param b the B input geometry
+     * @param pred the topological predicate
+     * @return true if the topological relationship is satisfied
+     */
+    static bool relate(const Geometry* a, const Geometry* b, TopologyPredicate& pred);
+
+    /**
+     * Tests whether the topological relationship between two geometries
+     * satisfies a topological predicate,
+     * using a given BoundaryNodeRule.
+     *
+     * @param a the A input geometry
+     * @param b the B input geometry
+     * @param pred the topological predicate
+     * @param bnRule the Boundary Node Rule to use
+     * @return true if the topological relationship is satisfied
+     */
+    static bool relate(const Geometry* a, const Geometry* b, TopologyPredicate& pred, const BoundaryNodeRule& bnRule);
+
+    /**
+     * Tests whether the topological relationship to a geometry
+     * matches a DE-9IM matrix pattern.
+     *
+     * @param a the A input geometry
+     * @param b the B input geometry
+     * @param imPattern the DE-9IM pattern to match
+     * @return true if the geometries relationship matches the DE-9IM pattern
+     *
+     * @see IntersectionMatrixPattern
+     */
+    static bool relate(const Geometry* a, const Geometry* b, const std::string& imPattern);
+
+    /**
+     * Computes the DE-9IM matrix
+     * for the topological relationship between two geometries.
+     *
+     * @param a the A input geometry
+     * @param b the B input geometry
+     * @return the DE-9IM matrix for the topological relationship
+     */
+    static IntersectionMatrix relate(const Geometry* a, const Geometry* b);
+
+    /**
+     * Computes the DE-9IM matrix
+     * for the topological relationship between two geometries.
+     *
+     * @param a the A input geometry
+     * @param b the B input geometry
+     * @param bnRule the Boundary Node Rule to use
+     * @return the DE-9IM matrix for the relationship
+     */
+    static IntersectionMatrix relate(const Geometry* a, const Geometry* b, const BoundaryNodeRule& bnRule);
+
+    /**
+     * Creates a prepared RelateNG instance to optimize the
+     * evaluation of relationships against a single geometry.
+     *
+     * @param a the A input geometry
+     * @return a prepared instance
+     */
+    static std::unique_ptr<RelateNG> prepare(const Geometry* a);
+
+    /**
+     * Creates a prepared RelateNG instance to optimize the
+     * computation of predicates against a single geometry,
+     * using a given BoundaryNodeRule.
+     *
+     * @param a the A input geometry
+     * @param bnRule the required BoundaryNodeRule
+     * @return a prepared instance
+     */
+    static std::unique_ptr<RelateNG> prepare(const Geometry* a, const BoundaryNodeRule& bnRule);
+
+
+    /**
+     * Computes the DE-9IM matrix for the topological relationship to a geometry.
+     *
+     * @param b the B geometry to test against
+     * @return the DE-9IM matrix
+     */
+    IntersectionMatrix evaluate(const Geometry* b);
+
+
+    /**
+     * Tests whether the topological relationship to a geometry
+     * matches a DE-9IM matrix pattern.
+     *
+     * @param b the B geometry to test against
+     * @param imPattern the DE-9IM pattern to match
+     * @return true if the geometries' topological relationship matches the DE-9IM pattern
+     *
+     * @see IntersectionMatrixPattern
+     */
+    bool evaluate(const Geometry* b, const std::string& imPattern);
+
+    /**
+     * Tests whether the topological relationship to a geometry
+     * satisfies a topology predicate.
+     *
+     * @param b the B geometry to test against
+     * @param predicate the topological predicate
+     * @return true if the predicate is satisfied
+     */
+    bool evaluate(const Geometry* b, TopologyPredicate& predicate);
+
+
+
+};
+
+} // namespace geos.operation.relateng
+} // namespace geos.operation
+} // namespace geos
+
diff --git a/include/geos/operation/relateng/RelateSegmentString.h b/include/geos/operation/relateng/RelateSegmentString.h
index d16ddb35b..556ac4d37 100644
--- a/include/geos/operation/relateng/RelateSegmentString.h
+++ b/include/geos/operation/relateng/RelateSegmentString.h
@@ -63,8 +63,6 @@ private:
     const RelateGeometry* m_inputGeom;
     const Geometry* m_parentPolygonal = nullptr;
 
-    std::unique_ptr<CoordinateSequence> csStore;
-
     // Constructor
     RelateSegmentString(
         const CoordinateSequence* pts,
@@ -73,8 +71,7 @@ private:
         int id,
         int ringId,
         const Geometry* poly,
-        const RelateGeometry* inputGeom,
-        bool orient)
+        const RelateGeometry* inputGeom)
         : BasicSegmentString(const_cast<CoordinateSequence*>(pts), nullptr)
         , m_isA(isA)
         , m_dimension(dimension)
@@ -82,21 +79,15 @@ private:
         , m_ringId(ringId)
         , m_inputGeom(inputGeom)
         , m_parentPolygonal(poly)
-        {
-            bool requireCW = (ringId == 0);
-            if (orient)
-                orientAndRemoveRepeated(requireCW);
-            else
-                removeRepeated();
-        }
+        {}
 
 
     // Methods
 
-    static std::unique_ptr<RelateSegmentString> createSegmentString(
+    static const RelateSegmentString* createSegmentString(
         const CoordinateSequence* pts,
         bool isA, int dim, int elementId, int ringId,
-        const Geometry* poly, const RelateGeometry* parent, bool orient);
+        const Geometry* poly, const RelateGeometry* parent);
 
     /**
      *
@@ -119,24 +110,18 @@ private:
         std::size_t segIndex,
         const CoordinateXY* pt) const;
 
-    void orientAndRemoveRepeated(bool orientCW);
-
-    void removeRepeated();
-
 
 public:
 
-    static std::unique_ptr<RelateSegmentString>
-    createLine(
+    static const RelateSegmentString* createLine(
         const CoordinateSequence* pts,
         bool isA, int elementId,
-        const RelateGeometry* parent, bool orient = false);
+        const RelateGeometry* parent);
 
-    static std::unique_ptr<RelateSegmentString>
-    createRing(
+    static const RelateSegmentString* createRing(
         const CoordinateSequence* pts,
         bool isA, int elementId, int ringId,
-        const Geometry* poly, const RelateGeometry* parent, bool orient = true);
+        const Geometry* poly, const RelateGeometry* parent);
 
     bool isA() const;
 
diff --git a/include/geos/operation/relateng/TopologyComputer.h b/include/geos/operation/relateng/TopologyComputer.h
index f4709940b..cad2a1bc3 100644
--- a/include/geos/operation/relateng/TopologyComputer.h
+++ b/include/geos/operation/relateng/TopologyComputer.h
@@ -52,9 +52,9 @@ class GEOS_DLL TopologyComputer {
 private:
 
     // Members
-    TopologyPredicate* predicate;
-    RelateGeometry* geomA;
-    RelateGeometry* geomB;
+    TopologyPredicate& predicate;
+    RelateGeometry& geomA;
+    RelateGeometry& geomB;
     std::map<CoordinateXY, NodeSections*> nodeMap;
     std::deque<std::unique_ptr<NodeSections>> nodeSectionsStore;
 
@@ -67,7 +67,7 @@ private:
 
     void initExteriorEmpty(bool geomNonEmpty);
 
-    RelateGeometry* getGeometry(bool isA) const;
+    RelateGeometry& getGeometry(bool isA) const;
 
     void updateDim(Location locA, Location locB, int dimension);
 
@@ -104,11 +104,11 @@ private:
 
     void addNodeSections(NodeSection* ns0, NodeSection* ns1);
 
-    void addLineEndOnPoint(bool isLineA, Location locLineEnd, Location locPoint);
+    void addLineEndOnPoint(bool isLineA, Location locLineEnd, Location locPoint, const CoordinateXY* pt);
 
-    void addLineEndOnLine(bool isLineA, Location locLineEnd, Location locLine);
+    void addLineEndOnLine(bool isLineA, Location locLineEnd, Location locLine, const CoordinateXY* pt);
 
-    void addLineEndOnArea(bool isLineA, Location locLineEnd, Location locArea);
+    void addLineEndOnArea(bool isLineA, Location locLineEnd, Location locArea, const CoordinateXY* pt);
 
     /**
      * Updates topology for an area vertex (in Interior or on Boundary)
@@ -121,9 +121,9 @@ private:
      * @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 addAreaVertexOnPoint(bool isAreaA, Location locArea, const CoordinateXY* pt);
 
-    void addAreaVertexOnLine(bool isAreaA, Location locArea, Location locTarget);
+    void addAreaVertexOnLine(bool isAreaA, Location locArea, Location locTarget, const CoordinateXY* pt);
 
     void evaluateNode(NodeSections* nodeSections);
 
@@ -136,9 +136,9 @@ private:
 public:
 
         TopologyComputer(
-            TopologyPredicate* p_predicate,
-            RelateGeometry* p_geomA,
-            RelateGeometry* p_geomB)
+            TopologyPredicate& p_predicate,
+            RelateGeometry& p_geomA,
+            RelateGeometry& p_geomB)
             : predicate(p_predicate)
             , geomA(p_geomA)
             , geomB(p_geomB)
@@ -180,13 +180,13 @@ public:
 
     void addIntersection(NodeSection* a, NodeSection* b);
 
-    void addPointOnPointInterior();
+    void addPointOnPointInterior(const CoordinateXY* pt);
 
-    void addPointOnPointExterior(bool isGeomA);
+    void addPointOnPointExterior(bool isGeomA, const CoordinateXY* pt);
 
-    void addPointOnGeometry(bool isA, Location locTarget, int dimTarget);
+    void addPointOnGeometry(bool isA, Location locTarget, int dimTarget, const CoordinateXY* pt);
 
-    void addLineEndOnGeometry(bool isLineA, Location locLineEnd, Location locTarget, int dimTarget);
+    void addLineEndOnGeometry(bool isLineA, Location locLineEnd, Location locTarget, int dimTarget, const CoordinateXY* pt);
 
     /**
      * Adds topology for an area vertex interaction with a target geometry element.
@@ -204,9 +204,9 @@ public:
      * @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 addAreaVertex(bool isAreaA, Location locArea, Location locTarget, int dimTarget, const CoordinateXY* pt);
 
-    void addAreaVertexOnArea(bool isAreaA, Location locArea, Location locTarget);
+    void addAreaVertexOnArea(bool isAreaA, Location locArea, Location locTarget, const CoordinateXY* pt);
 
     void evaluateNodes();
 
diff --git a/include/geos/operation/relateng/TopologyPredicate.h b/include/geos/operation/relateng/TopologyPredicate.h
index ce6e27e47..6b594dbca 100644
--- a/include/geos/operation/relateng/TopologyPredicate.h
+++ b/include/geos/operation/relateng/TopologyPredicate.h
@@ -191,6 +191,13 @@ public:
     virtual void updateDimension(Location locA, Location locB, int dimension) = 0;
 
 
+    friend std::ostream&
+    operator<<(std::ostream& os, const TopologyPredicate& ns)
+    {
+        os << ns.name();
+        return os;
+    }
+
 };
 
 } // namespace geos.operation.relateng
diff --git a/src/geom/GeometryCollection.cpp b/src/geom/GeometryCollection.cpp
index b1497eb4b..80ecbccae 100644
--- a/src/geom/GeometryCollection.cpp
+++ b/src/geom/GeometryCollection.cpp
@@ -456,23 +456,6 @@ GeometryCollection::reverseImpl() const
     return getFactory()->createGeometryCollection(std::move(reversed)).release();
 }
 
-void
-GeometryCollection::getAllGeometries(std::vector<const Geometry*>& geoms) const
-{
-    for (std::size_t i = 0; i < getNumGeometries(); i++) {
-        const Geometry* subGeom = getGeometryN(i);
-        if (subGeom == nullptr)
-            continue;
-
-        if (subGeom->isCollection()) {
-            const GeometryCollection* subColl = static_cast<const GeometryCollection*>(subGeom);
-            subColl->getAllGeometries(geoms);
-            continue;
-        }
-
-        geoms.push_back(subGeom);
-    }
-}
 
 
 } // namespace geos::geom
diff --git a/src/operation/relateng/EdgeSetIntersector.cpp b/src/operation/relateng/EdgeSetIntersector.cpp
index 826ee94ca..37b3600b8 100644
--- a/src/operation/relateng/EdgeSetIntersector.cpp
+++ b/src/operation/relateng/EdgeSetIntersector.cpp
@@ -38,19 +38,19 @@ namespace relateng {  // geos.operation.relateng
 
 /* private */
 void
-EdgeSetIntersector::addEdges(std::vector<std::unique_ptr<RelateSegmentString>>& segStrings)
+EdgeSetIntersector::addEdges(std::vector<const SegmentString*>& segStrings)
 {
-    for (auto& ss : segStrings) {
+    for (const SegmentString* ss : segStrings) {
         addToIndex(ss);
     }
 }
 
 /* private */
 void
-EdgeSetIntersector::addToIndex(std::unique_ptr<RelateSegmentString>& segStr)
+EdgeSetIntersector::addToIndex(const SegmentString* segStr)
 {
     std::vector<MonotoneChain> segChains;
-    MonotoneChainBuilder::getChains(segStr->getCoordinates(), segStr.get(), segChains);
+    MonotoneChainBuilder::getChains(segStr->getCoordinates(), const_cast<SegmentString*>(segStr), segChains);
 
     for (MonotoneChain& mc : segChains) {
         if (envelope == nullptr || envelope->intersects(mc.getEnvelope())) {
diff --git a/src/operation/relateng/RelateGeometry.cpp b/src/operation/relateng/RelateGeometry.cpp
index d31acf30e..7fd7998fa 100644
--- a/src/operation/relateng/RelateGeometry.cpp
+++ b/src/operation/relateng/RelateGeometry.cpp
@@ -13,7 +13,7 @@
  *
  **********************************************************************/
 
-#include <geos/algorithm/BoundaryNodeRule.h>
+#include <geos/algorithm/Orientation.h>
 #include <geos/geom/Coordinate.h>
 #include <geos/geom/Dimension.h>
 #include <geos/geom/Envelope.h>
@@ -27,16 +27,22 @@
 #include <geos/geom/Polygon.h>
 #include <geos/geom/util/ComponentCoordinateExtracter.h>
 #include <geos/geom/util/PointExtracter.h>
+#include <geos/geom/util/GeometryLister.h>
 #include <geos/operation/relateng/RelateGeometry.h>
 #include <geos/operation/relateng/RelateSegmentString.h>
 #include <geos/operation/relateng/DimensionLocation.h>
-
+#include <geos/operation/valid/RepeatedPointRemover.h>
 
 #include <sstream>
 
 
 using geos::algorithm::BoundaryNodeRule;
+using geos::algorithm::Orientation;
 using namespace geos::geom;
+using geos::geom::util::ComponentCoordinateExtracter;
+using geos::geom::util::GeometryLister;
+using geos::geom::util::PointExtracter;
+using geos::operation::valid::RepeatedPointRemover;
 
 
 namespace geos {      // geos
@@ -90,8 +96,7 @@ RelateGeometry::analyzeDimensions()
     }
     //-- analyze a (possibly mixed type) collection
     std::vector<const Geometry*> elems;
-    const GeometryCollection* col = static_cast<const GeometryCollection*>(geom);
-    col->getAllGeometries(elems);
+    GeometryLister::list(geom, elems);
     for (const Geometry* elem : elems)
     {
         if (elem->isEmpty())
@@ -117,8 +122,7 @@ bool
 RelateGeometry::isZeroLength(const Geometry* geom)
 {
     std::vector<const Geometry*> elems;
-    const GeometryCollection* col = static_cast<const GeometryCollection*>(geom);
-    col->getAllGeometries(elems);
+    GeometryLister::list(geom, elems);
     for (const Geometry* elem : elems) {
         if (elem->getGeometryTypeId() == GEOS_LINESTRING) {
             if (! isZeroLength(static_cast<const LineString*>(elem))) {
@@ -323,7 +327,7 @@ RelateGeometry::createUniquePoints()
 {
     //-- only called on P geometries
     std::vector<const CoordinateXY*> pts;
-    geom::util::ComponentCoordinateExtracter::getCoordinates(*geom, pts);
+    ComponentCoordinateExtracter::getCoordinates(*geom, pts);
     Coordinate::ConstXYSet set(pts.begin(), pts.end());
     return set;
 }
@@ -352,10 +356,11 @@ RelateGeometry::getEffectivePoints()
 
 
 /* public */
-std::vector<std::unique_ptr<RelateSegmentString>>
+std::vector<const SegmentString*>
 RelateGeometry::extractSegmentStrings(bool isA, const Envelope* env)
 {
-    std::vector<std::unique_ptr<RelateSegmentString>> segStrings;
+    std::vector<const SegmentString*> segStrings;
+    segStringStore.clear();
     extractSegmentStrings(isA, env, geom, segStrings);
     return segStrings;
 }
@@ -365,7 +370,7 @@ RelateGeometry::extractSegmentStrings(bool isA, const Envelope* env)
 void
 RelateGeometry::extractSegmentStrings(bool isA,
     const Envelope* env, const Geometry* p_geom,
-    std::vector<std::unique_ptr<RelateSegmentString>>& segStrings)
+    std::vector<const SegmentString*>& segStrings)
 {
     //-- record if parent is MultiPolygon
     const MultiPolygon* parentPolygonal = nullptr;
@@ -390,7 +395,7 @@ void
 RelateGeometry::extractSegmentStringsFromAtomic(bool isA,
     const Geometry* p_geom, const MultiPolygon* parentPolygonal,
     const Envelope* env,
-    std::vector<std::unique_ptr<RelateSegmentString>>& segStrings)
+    std::vector<const SegmentString*>& segStrings)
 {
     if (p_geom->isEmpty())
         return;
@@ -402,9 +407,14 @@ RelateGeometry::extractSegmentStringsFromAtomic(bool isA,
     elementId++;
     if (p_geom->getGeometryTypeId() == GEOS_LINESTRING) {
         const LineString* line = static_cast<const LineString*>(p_geom);
-        auto cs = line->getCoordinatesRO();
+        /*
+         * Condition the input Coordinate sequence so that it has no repeated points.
+         * This requires taking a copy which removeRepeated does behind the scenes and stores in csStore.
+         */
+        const CoordinateSequence* cs = removeRepeated(line->getCoordinatesRO());
         auto ss = RelateSegmentString::createLine(cs, isA, elementId, this);
-        segStrings.emplace_back(ss.release());
+        segStringStore.emplace_back(ss);
+        segStrings.push_back(ss);
     }
     else if (p_geom->getGeometryTypeId() == GEOS_POLYGON) {
         const Polygon* poly = static_cast<const Polygon*>(p_geom);
@@ -426,16 +436,23 @@ void
 RelateGeometry::extractRingToSegmentString(bool isA,
     const LinearRing* ring, int ringId, const Envelope* env,
     const Geometry* parentPoly,
-    std::vector<std::unique_ptr<RelateSegmentString>>& segStrings)
+    std::vector<const SegmentString*>& segStrings)
 {
     if (ring->isEmpty())
         return;
     if (env != nullptr && ! env->intersects(ring->getEnvelopeInternal()))
         return;
 
-    const CoordinateSequence* pts = ring->getCoordinatesRO();
-    auto ss = RelateSegmentString::createRing(pts, isA, elementId, ringId, parentPoly, this);
-    segStrings.emplace_back(ss.release());
+    /*
+     * Condition the input Coordinate sequence so that it has no repeated points
+     * and is oriented in a deterministic way. This requires taking a copy which
+     * orientAndRemoveRepeated does behind the scenes and stores in csStore.
+     */
+    bool requireCW = (ringId == 0);
+    const CoordinateSequence* cs = orientAndRemoveRepeated(ring->getCoordinatesRO(), requireCW);
+    auto ss = RelateSegmentString::createRing(cs, isA, elementId, ringId, parentPoly, this);
+    segStringStore.emplace_back(ss);
+    segStrings.push_back(ss);
 }
 
 
@@ -457,6 +474,51 @@ operator<<(std::ostream& os, const RelateGeometry& rg)
 
 
 
+/* private */
+const CoordinateSequence *
+RelateGeometry::orientAndRemoveRepeated(const CoordinateSequence *seq, bool orientCW)
+{
+    bool isFlipped = (orientCW == Orientation::isCCW(seq));
+    bool hasRepeated = seq->hasRepeatedPoints();
+    /* Already conditioned */
+    if (!isFlipped && !hasRepeated) {
+        return seq;
+    }
+
+    if (hasRepeated) {
+        auto deduped = RepeatedPointRemover::removeRepeatedPoints(seq);
+        if (isFlipped)
+            deduped->reverse();
+        CoordinateSequence* cs = deduped.release();
+        csStore.emplace_back(cs);
+        return cs;
+    }
+
+    if (isFlipped) {
+        auto reversed = seq->clone();
+        reversed->reverse();
+        CoordinateSequence* cs = reversed.release();
+        csStore.emplace_back(cs);
+        return cs;
+    }
+
+    return seq;
+}
+
+/* private */
+const CoordinateSequence *
+RelateGeometry::removeRepeated(const CoordinateSequence *seq)
+{
+    bool hasRepeated = seq->hasRepeatedPoints();
+    if (!hasRepeated)
+        return seq;
+    auto deduped = RepeatedPointRemover::removeRepeatedPoints(seq);
+    CoordinateSequence* cs = deduped.release();
+    csStore.emplace_back(cs);
+    return cs;
+}
+
+
 } // namespace geos.operation.overlayng
 } // namespace geos.operation
 } // namespace geos
diff --git a/src/operation/relateng/RelateNG.cpp b/src/operation/relateng/RelateNG.cpp
new file mode 100644
index 000000000..6b7a83f67
--- /dev/null
+++ b/src/operation/relateng/RelateNG.cpp
@@ -0,0 +1,523 @@
+/**********************************************************************
+ *
+ * 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/BoundaryNodeRule.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/Dimension.h>
+#include <geos/geom/Envelope.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/IntersectionMatrix.h>
+#include <geos/geom/LineString.h>
+#include <geos/geom/LinearRing.h>
+#include <geos/geom/Location.h>
+#include <geos/geom/Point.h>
+#include <geos/geom/Polygon.h>
+#include <geos/geom/prep/PreparedGeometry.h>
+#include <geos/geom/util/GeometryLister.h>
+#include <geos/noding/MCIndexSegmentSetMutualIntersector.h>
+#include <geos/noding/SegmentString.h>
+#include <geos/operation/relateng/DimensionLocation.h>
+#include <geos/operation/relateng/EdgeSegmentIntersector.h>
+#include <geos/operation/relateng/EdgeSetIntersector.h>
+#include <geos/operation/relateng/RelateGeometry.h>
+#include <geos/operation/relateng/RelateMatrixPredicate.h>
+#include <geos/operation/relateng/RelateNG.h>
+#include <geos/operation/relateng/RelatePredicate.h>
+#include <geos/operation/relateng/RelateSegmentString.h>
+#include <geos/operation/relateng/TopologyComputer.h>
+#include <geos/operation/relateng/TopologyPredicate.h>
+
+#include <sstream>
+
+
+using namespace geos::geom;
+using geos::algorithm::BoundaryNodeRule;
+using geos::noding::MCIndexSegmentSetMutualIntersector;
+using geos::noding::SegmentString;
+using geos::geom::util::GeometryLister;
+
+namespace geos {      // geos
+namespace operation { // geos.operation
+namespace relateng {  // geos.operation.relateng
+
+
+#define GEOM_A RelateGeometry::GEOM_A
+#define GEOM_B RelateGeometry::GEOM_B
+
+
+/* public static */
+bool
+RelateNG::relate(const Geometry* a, const Geometry* b, TopologyPredicate& pred)
+{
+    RelateNG rng(a, false);
+    return rng.evaluate(b, pred);
+}
+
+
+/* public static */
+bool
+RelateNG::relate(const Geometry* a, const Geometry* b, TopologyPredicate& pred, const BoundaryNodeRule& bnRule)
+{
+    RelateNG rng(a, false, bnRule);
+    return rng.evaluate(b, pred);
+}
+
+
+/* public static */
+bool
+RelateNG::relate(const Geometry* a, const Geometry* b, const std::string& imPattern)
+{
+    RelateNG rng(a, false);
+    return rng.evaluate(b, imPattern);
+}
+
+
+/* public static */
+IntersectionMatrix
+RelateNG::relate(const Geometry* a, const Geometry* b)
+{
+    RelateNG rng(a, false);
+    return rng.evaluate(b);
+}
+
+
+/* public static */
+IntersectionMatrix
+RelateNG::relate(const Geometry* a, const Geometry* b, const BoundaryNodeRule& bnRule)
+{
+    RelateNG rng(a, false, bnRule);
+    return rng.evaluate(b);
+}
+
+
+/* public static */
+std::unique_ptr<RelateNG>
+RelateNG::prepare(const Geometry* a)
+{
+    return std::unique_ptr<RelateNG>(new RelateNG(a, true));
+}
+
+
+/* public static */
+std::unique_ptr<RelateNG>
+RelateNG::prepare(const Geometry* a, const BoundaryNodeRule& bnRule)
+{
+    return std::unique_ptr<RelateNG>(new RelateNG(a, true, bnRule));
+}
+
+
+/* public */
+IntersectionMatrix
+RelateNG::evaluate(const Geometry* b)
+{
+    RelateMatrixPredicate rel;
+    evaluate(b, rel);
+    return rel.getIM();
+}
+
+
+/* public */
+bool
+RelateNG::evaluate(const Geometry* b, const std::string& imPattern)
+{
+    auto predicate = RelatePredicate::matches(imPattern);
+    return evaluate(b, *predicate);
+}
+
+
+/* public */
+bool
+RelateNG::evaluate(const Geometry* b, TopologyPredicate& predicate)
+{
+    //-- fast envelope checks
+    if (! hasRequiredEnvelopeInteraction(b, predicate)) {
+        return false;
+    }
+
+    RelateGeometry geomB(b, boundaryNodeRule);
+
+    if (geomA.isEmpty() && geomB.isEmpty()) {
+        //TODO: what if predicate is disjoint?  Perhaps use result on disjoint envs?
+        return finishValue(predicate);
+    }
+    int dimA = geomA.getDimensionReal();
+    int dimB = geomB.getDimensionReal();
+
+    //-- check if predicate is determined by dimension or envelope
+    predicate.init(dimA, dimB);
+    if (predicate.isKnown())
+        return finishValue(predicate);
+
+    predicate.init(*(geomA.getEnvelope()), *(geomB.getEnvelope()));
+    if (predicate.isKnown())
+        return finishValue(predicate);
+
+    TopologyComputer topoComputer(predicate, geomA, geomB);
+
+    //-- optimized P/P evaluation
+    if (dimA == Dimension::P && dimB == Dimension::P) {
+        computePP(geomB, topoComputer);
+        topoComputer.finish();
+        return topoComputer.getResult();
+    }
+
+    //-- test points against (potentially) indexed geometry first
+    computeAtPoints(geomB, GEOM_B, geomA, topoComputer);
+    if (topoComputer.isResultKnown()) {
+        return topoComputer.getResult();
+    }
+    computeAtPoints(geomA, GEOM_A, geomB, topoComputer);
+    if (topoComputer.isResultKnown()) {
+        return topoComputer.getResult();
+    }
+
+    if (geomA.hasEdges() && geomB.hasEdges()) {
+        computeAtEdges(geomB, topoComputer);
+    }
+
+    //-- after all processing, set remaining unknown values in IM
+    topoComputer.finish();
+    return topoComputer.getResult();
+}
+
+
+/* private */
+bool
+RelateNG::hasRequiredEnvelopeInteraction(const Geometry* b, TopologyPredicate& predicate)
+{
+    const Envelope* envB = b->getEnvelopeInternal();
+    bool isInteracts = false;
+    if (predicate.requireCovers(GEOM_A)) {
+        if (! geomA.getEnvelope()->covers(envB)) {
+            return false;
+        }
+        isInteracts = true;
+    }
+    else if (predicate.requireCovers(GEOM_B)) {
+        if (! envB->covers(geomA.getEnvelope())) {
+            return false;
+        }
+        isInteracts = true;
+    }
+    if (! isInteracts
+        && predicate.requireInteraction()
+        && ! geomA.getEnvelope()->intersects(envB)) {
+        return false;
+    }
+    return true;
+}
+
+/* private */
+bool
+RelateNG::finishValue(TopologyPredicate& predicate)
+{
+    predicate.finish();
+    return predicate.value();
+}
+
+
+/* private */
+void
+RelateNG::computePP(RelateGeometry& geomB, TopologyComputer& topoComputer)
+{
+    Coordinate::ConstXYSet& ptsA = geomA.getUniquePoints();
+    //TODO: only query points in interaction extent?
+    Coordinate::ConstXYSet& ptsB = geomB.getUniquePoints();
+
+    uint32_t numBinA = 0;
+    for (const CoordinateXY* ptB : ptsB) {
+        auto it = ptsA.find(ptB);
+        bool found = (it != ptsA.end());
+        if (found) {
+            numBinA++;
+            topoComputer.addPointOnPointInterior(ptB);
+        }
+        else {
+            topoComputer.addPointOnPointExterior(GEOM_B, ptB);
+        }
+        if (topoComputer.isResultKnown()) {
+            return;
+        }
+    }
+    /**
+     * If number of matched B points is less than size of A,
+     * there must be at least one A point in the exterior of B
+     */
+    if (numBinA < ptsA.size()) {
+        //TODO: determine actual exterior point?
+        topoComputer.addPointOnPointExterior(GEOM_A, nullptr);
+    }
+}
+
+
+/* private */
+void
+RelateNG::computeAtPoints(
+    RelateGeometry& geom, bool isA,
+    RelateGeometry& geomTarget, TopologyComputer& topoComputer)
+{
+    bool isResultKnown = false;
+    isResultKnown = computePoints(geom, isA, geomTarget, topoComputer);
+    if (isResultKnown)
+        return;
+
+    /**
+     * Performance optimization: only check points against target
+     * if it has areas OR if the predicate requires checking for
+     * exterior interaction.
+     * In particular, this avoids testing line ends against lines
+     * for the intersects predicate (since these are checked
+     * during segment/segment intersection checking anyway).
+     * Checking points against areas is necessary, since the input
+     * linework is disjoint if one input lies wholly inside an area,
+     * so segment intersection checking is not sufficient.
+     */
+    bool checkDisjointPoints = geomTarget.hasDimension(Dimension::A)
+                                || topoComputer.isExteriorCheckRequired(isA);
+    if (! checkDisjointPoints)
+        return;
+
+    isResultKnown = computeLineEnds(geom, isA, geomTarget, topoComputer);
+    if (isResultKnown)
+        return;
+
+    computeAreaVertex(geom, isA, geomTarget, topoComputer);
+}
+
+
+/* private */
+bool
+RelateNG::computePoints(
+    RelateGeometry& geom, bool isA,
+    RelateGeometry& geomTarget, TopologyComputer& topoComputer)
+{
+    if (! geom.hasDimension(Dimension::P)) {
+        return false;
+    }
+
+    std::vector<const Point*> points = geom.getEffectivePoints();
+    for (const Point* point : points) {
+        //TODO: exit when all possible target locations (E,I,B) have been found?
+        if (point->isEmpty())
+            continue;
+
+        const CoordinateXY* pt = point->getCoordinate();
+        computePoint(isA, pt, geomTarget, topoComputer);
+        if (topoComputer.isResultKnown()) {
+            return true;
+        }
+    }
+    return false;
+}
+
+
+/* private */
+void
+RelateNG::computePoint(bool isA, const CoordinateXY* pt, RelateGeometry& geomTarget, TopologyComputer& topoComputer)
+{
+      int locDimTarget = geomTarget.locateWithDim(pt);
+      Location locTarget = DimensionLocation::location(locDimTarget);
+      int dimTarget = DimensionLocation::dimension(locDimTarget, topoComputer.getDimension(! isA));
+      topoComputer.addPointOnGeometry(isA, locTarget, dimTarget, pt);
+}
+
+
+/* private */
+bool
+RelateNG::computeLineEnds(
+    RelateGeometry& geom, bool isA,
+    RelateGeometry& geomTarget, TopologyComputer& topoComputer)
+{
+    if (! geom.hasDimension(Dimension::L)) {
+        return false;
+    }
+
+    bool hasExteriorIntersection = false;
+    std::vector<const Geometry*> elems;
+    GeometryLister::list(geom.getGeometry(), elems);
+    for (const Geometry* elem : elems) {
+        if (elem->isEmpty())
+            continue;
+
+        if (elem->getGeometryTypeId() == GEOS_LINESTRING) {
+            //-- once an intersection with target exterior is recorded, skip further known-exterior points
+            if (hasExteriorIntersection
+                && elem->getEnvelopeInternal()->disjoint(geomTarget.getEnvelope()))
+                continue;
+
+            const LineString* line = static_cast<const LineString*>(elem);
+            //TODO: add optimzation to skip disjoint elements once exterior point found
+            const CoordinateXY& e0 = line->getCoordinatesRO()->getAt(0);
+            hasExteriorIntersection |= computeLineEnd(geom, isA, &e0, geomTarget, topoComputer);
+            if (topoComputer.isResultKnown()) {
+                return true;
+            }
+
+            if (! line->isClosed()) {
+                const CoordinateXY& e1 = line->getCoordinatesRO()->getAt(line->getNumPoints() - 1);
+                hasExteriorIntersection |= computeLineEnd(geom, isA, &e1, geomTarget, topoComputer);
+                if (topoComputer.isResultKnown()) {
+                    return true;
+                }
+            }
+        //TODO: break when all possible locations have been found?
+        }
+    }
+    return false;
+}
+
+
+/* private */
+bool
+RelateNG::computeLineEnd(RelateGeometry& geom, bool isA, const CoordinateXY* pt,
+      RelateGeometry& geomTarget, TopologyComputer& topoComputer)
+{
+    Location locLineEnd = geom.locateLineEnd(pt);
+    int locDimTarget = geomTarget.locateWithDim(pt);
+    Location locTarget = DimensionLocation::location(locDimTarget);
+    int dimTarget = DimensionLocation::dimension(locDimTarget, topoComputer.getDimension(! isA));
+    topoComputer.addLineEndOnGeometry(isA, locLineEnd, locTarget, dimTarget, pt);
+    return locTarget == Location::EXTERIOR;
+}
+
+
+/* private */
+bool
+RelateNG::computeAreaVertex(RelateGeometry& geom, bool isA, RelateGeometry& geomTarget, TopologyComputer& topoComputer)
+{
+    if (! geom.hasDimension(Dimension::A)) {
+        return false;
+    }
+    //-- evaluate for line and area targets only, since points are handled in the reverse direction
+    if (geomTarget.getDimension() < Dimension::L)
+        return false;
+
+    bool hasExteriorIntersection = false;
+
+    std::vector<const Geometry*> elems;
+    GeometryLister::list(geom.getGeometry(), elems);
+    for (const Geometry* elem : elems) {
+        if (elem->isEmpty())
+            continue;
+
+        if (elem->getGeometryTypeId() == GEOS_POLYGON) {
+            //-- once an intersection with target exterior is recorded, skip further known-exterior points
+            if (hasExteriorIntersection && elem->getEnvelopeInternal()->disjoint(geomTarget.getEnvelope()))
+                continue;
+
+            const Polygon* poly = static_cast<const Polygon*>(elem);
+            hasExteriorIntersection |= computeAreaVertex(geom, isA, poly->getExteriorRing(), geomTarget, topoComputer);
+            if (topoComputer.isResultKnown()) {
+                return true;
+            }
+            for (std::size_t j = 0; j < poly->getNumInteriorRing(); j++) {
+                hasExteriorIntersection |= computeAreaVertex(geom, isA, poly->getInteriorRingN(j), geomTarget, topoComputer);
+                if (topoComputer.isResultKnown()) {
+                    return true;
+                }
+            }
+        }
+    }
+    return false;
+}
+
+
+/* private */
+bool
+RelateNG::computeAreaVertex(RelateGeometry& geom, bool isA, const LinearRing* ring, RelateGeometry& geomTarget, TopologyComputer& topoComputer)
+{
+    //TODO: use extremal (highest) point to ensure one is on boundary of polygon cluster
+    const CoordinateXY* pt = ring->getCoordinate();
+
+    Location locArea = geom.locateAreaVertex(pt);
+    int locDimTarget = geomTarget.locateWithDim(pt);
+    Location locTarget = DimensionLocation::location(locDimTarget);
+    int dimTarget = DimensionLocation::dimension(locDimTarget, topoComputer.getDimension(! isA));
+    topoComputer.addAreaVertex(isA, locArea, locTarget, dimTarget, pt);
+    return locTarget == Location::EXTERIOR;
+}
+
+
+/* private */
+void
+RelateNG::computeAtEdges(RelateGeometry& geomB, TopologyComputer& topoComputer)
+{
+    Envelope envInt;
+    geomA.getEnvelope()->intersection(*(geomB.getEnvelope()), envInt);
+    if (envInt.isNull())
+        return;
+
+    std::vector<const SegmentString*> edgesB = geomB.extractSegmentStrings(GEOM_B, &envInt);
+    EdgeSegmentIntersector intersector(topoComputer);
+
+    if (topoComputer.isSelfNodingRequired()) {
+        computeEdgesAll(edgesB, &envInt, intersector);
+    }
+    else {
+        computeEdgesMutual(edgesB, &envInt, intersector);
+    }
+    if (topoComputer.isResultKnown()) {
+        return;
+    }
+
+    topoComputer.evaluateNodes();
+}
+
+
+/* private */
+void
+RelateNG::computeEdgesAll(std::vector<const SegmentString*>& edgesB, const Envelope* envInt, EdgeSegmentIntersector& intersector)
+{
+    //TODO: find a way to reuse prepared index?
+    std::vector<const SegmentString*> edgesA = geomA.extractSegmentStrings(GEOM_A, envInt);
+
+    EdgeSetIntersector edgeInt(edgesA, edgesB, envInt);
+    edgeInt.process(intersector);
+}
+
+
+/* private */
+void
+RelateNG::computeEdgesMutual(std::vector<const SegmentString*>& edgesB, const Envelope* envInt, EdgeSegmentIntersector& intersector)
+{
+    //-- in prepared mode the A edge index is reused
+    if (edgeMutualInt == nullptr) {
+        const Envelope* envExtract = geomA.isPrepared() ? nullptr : envInt;
+        std::vector<const SegmentString*> edgesA = geomA.extractSegmentStrings(GEOM_A, envExtract);
+        edgeMutualInt.reset(new MCIndexSegmentSetMutualIntersector(envExtract));
+        edgeMutualInt->setBaseSegments(&edgesA);
+        edgeMutualInt->setSegmentIntersector(&intersector);
+
+    }
+
+    edgeMutualInt->process(&edgesB);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+} // namespace geos.operation.overlayng
+} // namespace geos.operation
+} // namespace geos
+
+
diff --git a/src/operation/relateng/RelateSegmentString.cpp b/src/operation/relateng/RelateSegmentString.cpp
index d6f0815a9..5dfb5e716 100644
--- a/src/operation/relateng/RelateSegmentString.cpp
+++ b/src/operation/relateng/RelateSegmentString.cpp
@@ -18,12 +18,9 @@
 #include <geos/geom/Coordinate.h>
 #include <geos/geom/CoordinateSequence.h>
 #include <geos/geom/Geometry.h>
-
-// #include <geos/io/WKTWriter.h>
 #include <geos/operation/relateng/NodeSection.h>
 #include <geos/operation/relateng/RelateGeometry.h>
 #include <geos/operation/relateng/RelateSegmentString.h>
-#include <geos/operation/valid/RepeatedPointRemover.h>
 #include <sstream>
 
 
@@ -32,7 +29,6 @@ using geos::geom::CoordinateXY;
 using geos::geom::Dimension;
 using geos::geom::Geometry;
 using geos::algorithm::Orientation;
-using geos::operation::valid::RepeatedPointRemover;
 
 
 namespace geos {      // geos
@@ -41,35 +37,35 @@ namespace relateng {  // geos.operation.relateng
 
 
 /* public static */
-std::unique_ptr<RelateSegmentString>
+const RelateSegmentString*
 RelateSegmentString::createLine(
     const CoordinateSequence* pts,
     bool isA, int elementId,
-    const RelateGeometry* parent, bool orient)
+    const RelateGeometry* parent)
 {
-    return createSegmentString(pts, isA, Dimension::L, elementId, -1, nullptr, parent, orient);
+    return createSegmentString(pts, isA, Dimension::L, elementId, -1, nullptr, parent);
 }
 
 
 /* public static */
-std::unique_ptr<RelateSegmentString>
+const RelateSegmentString*
 RelateSegmentString::createRing(
     const CoordinateSequence* pts,
     bool isA, int elementId, int ringId,
-    const Geometry* poly, const RelateGeometry* parent, bool orient)
+    const Geometry* poly, const RelateGeometry* parent)
 {
-    return createSegmentString(pts, isA, Dimension::A, elementId, ringId, poly, parent, orient);
+    return createSegmentString(pts, isA, Dimension::A, elementId, ringId, poly, parent);
 }
 
 
 /* private static */
-std::unique_ptr<RelateSegmentString>
+const RelateSegmentString*
 RelateSegmentString::createSegmentString(
     const CoordinateSequence* pts,
     bool isA, int dim, int elementId, int ringId,
-    const Geometry* poly, const RelateGeometry* parent, bool orient)
+    const Geometry* poly, const RelateGeometry* parent)
 {
-    return std::unique_ptr<RelateSegmentString>(new RelateSegmentString(pts, isA, dim, elementId, ringId, poly, parent, orient));
+    return new RelateSegmentString(pts, isA, dim, elementId, ringId, poly, parent);
 }
 
 
@@ -180,41 +176,6 @@ RelateSegmentString::isContainingSegment(std::size_t segIndex, const CoordinateX
 }
 
 
-/* public */
-void
-RelateSegmentString::orientAndRemoveRepeated(bool orientCW)
-{
-    bool isFlipped = (orientCW == Orientation::isCCW(seq));
-    bool hasRepeated = seq->hasRepeatedPoints();
-    /* Already conditioned */
-    if (!isFlipped && !hasRepeated) {
-        return;
-    }
-
-    if (hasRepeated) {
-        csStore = RepeatedPointRemover::removeRepeatedPoints(seq);
-        if (isFlipped)
-            csStore->reverse();
-    }
-
-    if (isFlipped) {
-        csStore = seq->clone();
-        csStore->reverse();
-    }
-
-    seq = csStore.get();
-}
-
-/* public */
-void
-RelateSegmentString::removeRepeated()
-{
-    bool hasRepeated = seq->hasRepeatedPoints();
-    if (!hasRepeated)
-        return;
-    csStore = RepeatedPointRemover::removeRepeatedPoints(seq);
-    seq = csStore.get();
-}
 
 
 } // namespace geos.operation.overlayng
diff --git a/src/operation/relateng/TopologyComputer.cpp b/src/operation/relateng/TopologyComputer.cpp
index 37adebd6c..78dd304ea 100644
--- a/src/operation/relateng/TopologyComputer.cpp
+++ b/src/operation/relateng/TopologyComputer.cpp
@@ -47,8 +47,8 @@ namespace relateng {  // geos.operation.relateng
 void
 TopologyComputer::initExteriorDims()
 {
-    int dimRealA = geomA->getDimensionReal();
-    int dimRealB = geomB->getDimensionReal();
+    int dimRealA = geomA.getDimensionReal();
+    int dimRealB = geomB.getDimensionReal();
 
     /**
      * For P/L case, P exterior intersects L interior
@@ -98,7 +98,7 @@ TopologyComputer::initExteriorEmpty(bool geomNonEmpty)
             updateDim(geomNonEmpty, Location::INTERIOR, Location::EXTERIOR, Dimension::P);
             break;
         case Dimension::L:
-            if (getGeometry(geomNonEmpty)->hasBoundary()) {
+            if (getGeometry(geomNonEmpty).hasBoundary()) {
                 updateDim(geomNonEmpty, Location::BOUNDARY, Location::EXTERIOR, Dimension::P);
             }
             updateDim(geomNonEmpty, Location::INTERIOR, Location::EXTERIOR, Dimension::L);
@@ -112,7 +112,7 @@ TopologyComputer::initExteriorEmpty(bool geomNonEmpty)
 
 
 /* private */
-RelateGeometry*
+RelateGeometry&
 TopologyComputer::getGeometry(bool isA) const
 {
     return isA ? geomA : geomB;
@@ -123,7 +123,7 @@ TopologyComputer::getGeometry(bool isA) const
 int
 TopologyComputer::getDimension(bool isA) const
 {
-    return getGeometry(isA)->getDimension();
+    return getGeometry(isA).getDimension();
 }
 
 
@@ -141,9 +141,9 @@ 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();
+    if (geomA.isPointsOrPolygons()) return false;
+    if (geomB.isPointsOrPolygons()) return false;
+    return predicate.requireSelfNoding();
 }
 
 
@@ -151,7 +151,7 @@ TopologyComputer::isSelfNodingRequired() const
 bool
 TopologyComputer::isExteriorCheckRequired(bool isA) const
 {
-    return predicate->requireExteriorCheck(isA);
+    return predicate.requireExteriorCheck(isA);
 }
 
 
@@ -159,7 +159,7 @@ TopologyComputer::isExteriorCheckRequired(bool isA) const
 void
 TopologyComputer::updateDim(Location locA, Location locB, int dimension)
 {
-    predicate->updateDimension(locA, locB, dimension);
+    predicate.updateDimension(locA, locB, dimension);
 }
 
 
@@ -181,7 +181,7 @@ TopologyComputer::updateDim(bool isAB, Location loc1, Location loc2, int dimensi
 bool
 TopologyComputer::isResultKnown() const
 {
-    return predicate->isKnown();
+    return predicate.isKnown();
 }
 
 
@@ -189,7 +189,7 @@ TopologyComputer::isResultKnown() const
 bool
 TopologyComputer::getResult() const
 {
-    return predicate->value();
+    return predicate.value();
 }
 
 
@@ -197,7 +197,7 @@ TopologyComputer::getResult() const
 void
 TopologyComputer::finish()
 {
-    predicate->finish();
+    predicate.finish();
 }
 
 
@@ -261,8 +261,8 @@ 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());
+    Location locA = geomA.locateNode(&pt, a->getPolygonal());
+    Location locB = geomB.locateNode(&pt, b->getPolygonal());
     updateDim(locA, locB, Dimension::P);
 }
 
@@ -277,22 +277,25 @@ TopologyComputer::addNodeSections(NodeSection* ns0, NodeSection* ns1)
 
 /* public */
 void
-TopologyComputer::addPointOnPointInterior()
+TopologyComputer::addPointOnPointInterior(const CoordinateXY* pt)
 {
+    (void)pt;
     updateDim(Location::INTERIOR, Location::INTERIOR, Dimension::P);
 }
 
 /* public */
 void
-TopologyComputer::addPointOnPointExterior(bool isGeomA)
+TopologyComputer::addPointOnPointExterior(bool isGeomA, const CoordinateXY* pt)
 {
+    (void)pt;
     updateDim(isGeomA, Location::INTERIOR, Location::EXTERIOR, Dimension::P);
 }
 
 /* public */
 void
-TopologyComputer::addPointOnGeometry(bool isA, Location locTarget, int dimTarget)
+TopologyComputer::addPointOnGeometry(bool isA, Location locTarget, int dimTarget, const CoordinateXY* pt)
 {
+    (void)pt;
     updateDim(isA, Location::INTERIOR, locTarget, Dimension::P);
     switch (dimTarget) {
     case Dimension::P:
@@ -320,17 +323,18 @@ TopologyComputer::addPointOnGeometry(bool isA, Location locTarget, int dimTarget
 
 /* public */
 void
-TopologyComputer::addLineEndOnGeometry(bool isLineA, Location locLineEnd, Location locTarget, int dimTarget)
+TopologyComputer::addLineEndOnGeometry(bool isLineA, Location locLineEnd, Location locTarget, int dimTarget, const CoordinateXY* pt)
 {
+    (void)pt;
     switch (dimTarget) {
     case Dimension::P:
-        addLineEndOnPoint(isLineA, locLineEnd, locTarget);
+        addLineEndOnPoint(isLineA, locLineEnd, locTarget, pt);
         return;
     case Dimension::L:
-        addLineEndOnLine(isLineA, locLineEnd, locTarget);
+        addLineEndOnLine(isLineA, locLineEnd, locTarget, pt);
         return;
     case Dimension::A:
-        addLineEndOnArea(isLineA, locLineEnd, locTarget);
+        addLineEndOnArea(isLineA, locLineEnd, locTarget, pt);
         return;
     }
     throw IllegalStateException("Unknown target dimension: " + std::to_string(dimTarget));
@@ -338,15 +342,17 @@ TopologyComputer::addLineEndOnGeometry(bool isLineA, Location locLineEnd, Locati
 
 /* private */
 void
-TopologyComputer::addLineEndOnPoint(bool isLineA, Location locLineEnd, Location locPoint)
+TopologyComputer::addLineEndOnPoint(bool isLineA, Location locLineEnd, Location locPoint, const CoordinateXY* pt)
 {
+    (void)pt;
     updateDim(isLineA, locLineEnd, locPoint, Dimension::P);
 }
 
 /* private */
 void
-TopologyComputer::addLineEndOnLine(bool isLineA, Location locLineEnd, Location locLine)
+TopologyComputer::addLineEndOnLine(bool isLineA, Location locLineEnd, Location locLine, const CoordinateXY* pt)
 {
+    (void)pt;
     updateDim(isLineA, locLineEnd, locLine, Dimension::P);
     /**
      * When a line end is in the exterior, some length of the line interior
@@ -361,8 +367,9 @@ TopologyComputer::addLineEndOnLine(bool isLineA, Location locLineEnd, Location l
 
 /* private */
 void
-TopologyComputer::addLineEndOnArea(bool isLineA, Location locLineEnd, Location locArea)
+TopologyComputer::addLineEndOnArea(bool isLineA, Location locLineEnd, Location locArea, const CoordinateXY* pt)
 {
+    (void)pt;
     if (locArea == Location::BOUNDARY) {
         updateDim(isLineA, locLineEnd, locArea, Dimension::P);
     }
@@ -377,8 +384,9 @@ TopologyComputer::addLineEndOnArea(bool isLineA, Location locLineEnd, Location l
 
 /* public */
 void
-TopologyComputer::addAreaVertex(bool isAreaA, Location locArea, Location locTarget, int dimTarget)
+TopologyComputer::addAreaVertex(bool isAreaA, Location locArea, Location locTarget, int dimTarget, const CoordinateXY* pt)
 {
+    (void)pt;
     if (locTarget == Location::EXTERIOR) {
         updateDim(isAreaA, Location::INTERIOR, Location::EXTERIOR, Dimension::A);
         /**
@@ -397,13 +405,13 @@ TopologyComputer::addAreaVertex(bool isAreaA, Location locArea, Location locTarg
 
     switch (dimTarget) {
         case Dimension::P:
-            addAreaVertexOnPoint(isAreaA, locArea);
+            addAreaVertexOnPoint(isAreaA, locArea, pt);
             return;
         case Dimension::L:
-            addAreaVertexOnLine(isAreaA, locArea, locTarget);
+            addAreaVertexOnLine(isAreaA, locArea, locTarget, pt);
             return;
         case Dimension::A:
-            addAreaVertexOnArea(isAreaA, locArea, locTarget);
+            addAreaVertexOnArea(isAreaA, locArea, locTarget, pt);
             return;
     }
     throw IllegalStateException("Unknown target dimension: " + std::to_string(dimTarget));
@@ -412,8 +420,9 @@ TopologyComputer::addAreaVertex(bool isAreaA, Location locArea, Location locTarg
 
 /* private */
 void
-TopologyComputer::addAreaVertexOnPoint(bool isAreaA, Location locArea)
+TopologyComputer::addAreaVertexOnPoint(bool isAreaA, Location locArea, const CoordinateXY* pt)
 {
+    (void)pt;
     //-- Assert: locArea != EXTERIOR
     //-- Assert: locTarget == INTERIOR
     /**
@@ -437,8 +446,9 @@ TopologyComputer::addAreaVertexOnPoint(bool isAreaA, Location locArea)
 
 /* private */
 void
-TopologyComputer::addAreaVertexOnLine(bool isAreaA, Location locArea, Location locTarget)
+TopologyComputer::addAreaVertexOnLine(bool isAreaA, Location locArea, Location locTarget, const CoordinateXY* pt)
 {
+    (void)pt;
     //-- Assert: locArea != EXTERIOR
     /**
      * If an area vertex intersects a line, all we know is the
@@ -459,8 +469,9 @@ TopologyComputer::addAreaVertexOnLine(bool isAreaA, Location locArea, Location l
 
 /* public */
 void
-TopologyComputer::addAreaVertexOnArea(bool isAreaA, Location locArea, Location locTarget)
+TopologyComputer::addAreaVertexOnArea(bool isAreaA, Location locArea, Location locTarget, const CoordinateXY* pt)
 {
+    (void)pt;
     if (locTarget == Location::BOUNDARY) {
         if (locArea == Location::BOUNDARY) {
             //-- B/B topology is fully computed later by node analysis
@@ -513,8 +524,8 @@ 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));
+    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());
 }
diff --git a/tests/unit/geom/util/GeometryListerTest.cpp b/tests/unit/geom/util/GeometryListerTest.cpp
new file mode 100644
index 000000000..c63a3dad1
--- /dev/null
+++ b/tests/unit/geom/util/GeometryListerTest.cpp
@@ -0,0 +1,62 @@
+//
+// Test Suite for geos::geom::util::GeometryLister class.
+
+// tut
+#include <tut/tut.hpp>
+
+// geos
+#include <geos/geom/Geometry.h>
+#include <geos/io/WKTReader.h>
+#include <geos/geom/util/GeometryLister.h>
+
+// std
+#include <vector>
+
+namespace tut {
+//
+// Test Group
+//
+
+using namespace geos::geom::util;
+using namespace geos::geom;
+
+// Common data used by tests
+struct test_geometrylister_data {
+
+    geos::io::WKTReader _reader;
+
+    test_geometrylister_data() {};
+
+    void checkList(const std::string& wkt, std::size_t nSingletons)
+    {
+        std::unique_ptr<Geometry> geom(_reader.read(wkt));
+        std::vector<const Geometry*> elems;
+        GeometryLister::list(geom.get(), elems);
+        ensure_equals(elems.size(), nSingletons);
+    }
+};
+
+typedef test_group<test_geometrylister_data> group;
+typedef group::object object;
+
+group test_geometrylister_group("geos::geom::util::GeometryLister");
+
+//
+// Test Cases
+//
+
+template<>
+template<>
+void object::test<1> ()
+{
+    checkList("POINT(1 1)", 1);
+}
+
+template<>
+template<>
+void object::test<2> ()
+{
+    checkList("GEOMETRYCOLLECTION(MULTIPOINT(-117 33,-33 44),LINESTRING(0 0, 10 0),POINT(0 0),POLYGON((0 0, 10 0, 10 10, 0 10, 0 0)),GEOMETRYCOLLECTION(POINT(3 4)))", 6);
+}
+
+} // namespace tut
diff --git a/tests/unit/operation/relateng/RelateNGTest.cpp b/tests/unit/operation/relateng/RelateNGTest.cpp
new file mode 100644
index 000000000..4f45481db
--- /dev/null
+++ b/tests/unit/operation/relateng/RelateNGTest.cpp
@@ -0,0 +1,935 @@
+//
+// Test Suite for geos::operation::relateng::RelateNG class.
+
+#include <tut/tut.hpp>
+#include <utility.h>
+
+// geos
+#include <geos/io/WKTReader.h>
+// #include <geos/io/WKTWriter.h>
+#include <geos/geom/Geometry.h>
+#include <geos/operation/relateng/RelateNG.h>
+#include <geos/operation/relateng/RelatePredicate.h>
+#include <geos/operation/relateng/RelateMatrixPredicate.h>
+#include <geos/operation/relateng/IntersectionMatrixPattern.h>
+
+// std
+#include <memory>
+
+using namespace geos::geom;
+using namespace geos::operation::relateng;
+using geos::io::WKTReader;
+// using geos::io::WKTWriter;
+
+namespace tut {
+//
+// Test Group
+//
+
+// Common data used by all tests
+struct test_relateng_data {
+
+    WKTReader r;
+    // WKTWriter w;
+
+    void checkIntersectsDisjoint(const std::string& wkta, const std::string& wktb, bool expectedValue)
+    {
+        checkPredicate(*RelatePredicate::intersects(), wkta, wktb, expectedValue);
+        checkPredicate(*RelatePredicate::intersects(), wktb, wkta, expectedValue);
+        checkPredicate(*RelatePredicate::disjoint(), wkta, wktb, ! expectedValue);
+        checkPredicate(*RelatePredicate::disjoint(), wktb, wkta, ! expectedValue);
+    }
+
+    void checkContainsWithin(const std::string& wkta, const std::string& wktb, bool expectedValue)
+    {
+        checkPredicate(*RelatePredicate::contains(), wkta, wktb, expectedValue);
+        checkPredicate(*RelatePredicate::within(),   wktb, wkta, expectedValue);
+    }
+
+    void checkCoversCoveredBy(const std::string& wkta, const std::string& wktb, bool expectedValue)
+    {
+        checkPredicate(*RelatePredicate::covers(),    wkta, wktb, expectedValue);
+        checkPredicate(*RelatePredicate::coveredBy(), wktb, wkta, expectedValue);
+    }
+
+    void checkCrosses(const std::string& wkta, const std::string& wktb, bool expectedValue)
+    {
+        checkPredicate(*RelatePredicate::crosses(), wkta, wktb, expectedValue);
+        checkPredicate(*RelatePredicate::crosses(), wktb, wkta, expectedValue);
+    }
+
+    void checkOverlaps(const std::string& wkta, const std::string& wktb, bool expectedValue)
+    {
+        checkPredicate(*RelatePredicate::overlaps(), wkta, wktb, expectedValue);
+        checkPredicate(*RelatePredicate::overlaps(), wktb, wkta, expectedValue);
+    }
+
+    void checkTouches(const std::string& wkta, const std::string& wktb, bool expectedValue)
+    {
+        checkPredicate(*RelatePredicate::touches(), wkta, wktb, expectedValue);
+        checkPredicate(*RelatePredicate::touches(), wktb, wkta, expectedValue);
+    }
+
+    void checkEquals(const std::string& wkta, const std::string& wktb, bool expectedValue)
+    {
+        checkPredicate(*RelatePredicate::equalsTopo(), wkta, wktb, expectedValue);
+        checkPredicate(*RelatePredicate::equalsTopo(), wktb, wkta, expectedValue);
+    }
+
+    void checkRelate(const std::string& wkta, const std::string& wktb, const std::string expectedValue)
+    {
+        std::unique_ptr<Geometry> a = r.read(wkta);
+        std::unique_ptr<Geometry> b = r.read(wktb);
+        RelateMatrixPredicate pred;
+        // TopologyPredicate predTrace = trace(pred);
+        RelateNG::relate(a.get(), b.get(), pred);
+        std::string actualVal = pred.getIM().toString();
+        ensure_equals("checkRelate", actualVal, expectedValue);
+    }
+
+    void checkRelateMatches(const std::string& wkta, const std::string& wktb, const std::string pattern, bool expectedValue)
+    {
+        auto pred = RelatePredicate::matches(pattern);
+        checkPredicate(*pred, wkta, wktb, expectedValue);
+    }
+
+    void checkPredicate(TopologyPredicate& pred, const std::string& wkta, const std::string& wktb, bool expectedValue)
+    {
+        std::unique_ptr<Geometry> a = r.read(wkta);
+        std::unique_ptr<Geometry> b = r.read(wktb);
+        // TopologyPredicate predTrace = trace(pred);
+        bool actualVal = RelateNG::relate(a.get(), b.get(), pred);
+        if (actualVal != expectedValue) {
+            std::cerr << *a << " " << pred << " " << *b << " = " << actualVal << std::endl;
+        }
+        ensure_equals("checkPredicate", actualVal, expectedValue);
+    }
+
+};
+
+typedef test_group<test_relateng_data> group;
+typedef group::object object;
+
+group test_relateng_group("geos::operation::relateng::RelateNG");
+
+
+// testPointsDisjoint
+template<>
+template<>
+void object::test<1> ()
+{
+    std::string a = "POINT (0 0)";
+    std::string b = "POINT (1 1)";
+    checkIntersectsDisjoint(a, b, false);
+    checkContainsWithin(a, b, false);
+    checkEquals(a, b, false);
+    checkRelate(a, b, "FF0FFF0F2");
+}
+
+//======= P/P  =============
+
+// testPointsContained
+template<>
+template<>
+void object::test<2> ()
+{
+    std::string a = "MULTIPOINT (0 0, 1 1, 2 2)";
+    std::string b = "MULTIPOINT (1 1, 2 2)";
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, true);
+    checkEquals(a, b, false);
+    checkRelate(a, b, "0F0FFFFF2");
+}
+
+// testPointsEqual
+template<>
+template<>
+void object::test<3> ()
+{
+    std::string a = "MULTIPOINT (0 0, 1 1, 2 2)";
+    std::string b = "MULTIPOINT (0 0, 1 1, 2 2)";
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, true);
+    checkEquals(a, b, true);
+}
+
+// testValidateRelatePP_13
+template<>
+template<>
+void object::test<4> ()
+{
+    std::string a = "MULTIPOINT ((80 70), (140 120), (20 20), (200 170))";
+    std::string b = "MULTIPOINT ((80 70), (140 120), (80 170), (200 80))";
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, false);
+    checkContainsWithin(b, a, false);
+    checkCoversCoveredBy(a, b, false);
+    checkOverlaps(a, b, true);
+    checkTouches(a, b, false);
+}
+
+//======= L/P  =============
+
+// testLinePointContains
+template<>
+template<>
+void object::test<5> ()
+{
+    std::string a = "LINESTRING (0 0, 1 1, 2 2)";
+    std::string b = "MULTIPOINT (0 0, 1 1, 2 2)";
+    checkRelate(a, b, "0F10FFFF2");
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, true);
+    checkContainsWithin(b, a, false);
+    checkCoversCoveredBy(a, b, true);
+    checkCoversCoveredBy(b, a, false);
+}
+
+// testLinePointOverlaps
+template<>
+template<>
+void object::test<6> ()
+{
+    std::string a = "LINESTRING (0 0, 1 1)";
+    std::string b = "MULTIPOINT (0 0, 1 1, 2 2)";
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, false);
+    checkContainsWithin(b, a, false);
+    checkCoversCoveredBy(a, b, false);
+    checkCoversCoveredBy(b, a, false);
+}
+
+// testZeroLengthLinePoint
+template<>
+template<>
+void object::test<7> ()
+{
+    std::string a = "LINESTRING (0 0, 0 0)";
+    std::string b = "POINT (0 0)";
+    checkRelate(a, b, "0FFFFFFF2");
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, true);
+    checkContainsWithin(b, a, true);
+    checkCoversCoveredBy(a, b, true);
+    checkCoversCoveredBy(b, a, true);
+    checkEquals(a, b, true);
+}
+
+// testZeroLengthLineLine
+template<>
+template<>
+void object::test<8> ()
+{
+    std::string a = "LINESTRING (10 10, 10 10, 10 10)";
+    std::string b = "LINESTRING (10 10, 10 10)";
+    checkRelate(a, b, "0FFFFFFF2");
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, true);
+    checkContainsWithin(b, a, true);
+    checkCoversCoveredBy(a, b, true);
+    checkCoversCoveredBy(b, a, true);
+    checkEquals(a, b, true);
+}
+
+// tests bug involving checking for non-zero-length lines
+// testNonZeroLengthLinePoint
+template<>
+template<>
+void object::test<9> ()
+{
+    std::string a = "LINESTRING (0 0, 0 0, 9 9)";
+    std::string b = "POINT (1 1)";
+    checkRelate(a, b, "0F1FF0FF2");
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, true);
+    checkContainsWithin(b, a, false);
+    checkCoversCoveredBy(a, b, true);
+    checkCoversCoveredBy(b, a, false);
+    checkEquals(a, b, false);
+}
+
+// testLinePointIntAndExt
+template<>
+template<>
+void object::test<10> ()
+{
+    std::string a = "MULTIPOINT((60 60), (100 100))";
+    std::string b = "LINESTRING(40 40, 80 80)";
+    checkRelate(a, b, "0F0FFF102");
+}
+
+//======= L/L  =============
+
+// testLinesCrossProper
+template<>
+template<>
+void object::test<11> ()
+{
+    std::string a = "LINESTRING (0 0, 9 9)";
+    std::string b = "LINESTRING(0 9, 9 0)";
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, false);
+}
+
+// testLinesOverlap
+template<>
+template<>
+void object::test<12> ()
+{
+    std::string a = "LINESTRING (0 0, 5 5)";
+    std::string b = "LINESTRING(3 3, 9 9)";
+    checkIntersectsDisjoint(a, b, true);
+    checkTouches(a, b, false);
+    checkOverlaps(a, b, true);
+}
+
+// testLinesCrossVertex
+template<>
+template<>
+void object::test<13> ()
+{
+    std::string a = "LINESTRING (0 0, 8 8)";
+    std::string b = "LINESTRING(0 8, 4 4, 8 0)";
+    checkIntersectsDisjoint(a, b, true);
+}
+
+// testLinesTouchVertex
+template<>
+template<>
+void object::test<14> ()
+{
+    std::string a = "LINESTRING (0 0, 8 0)";
+    std::string b = "LINESTRING(0 8, 4 0, 8 8)";
+    checkIntersectsDisjoint(a, b, true);
+}
+
+// testLinesDisjointByEnvelope
+template<>
+template<>
+void object::test<15> ()
+{
+    std::string a = "LINESTRING (0 0, 9 9)";
+    std::string b = "LINESTRING(10 19, 19 10)";
+    checkIntersectsDisjoint(a, b, false);
+    checkContainsWithin(a, b, false);
+}
+
+// testLinesDisjoint
+template<>
+template<>
+void object::test<16> ()
+{
+    std::string a = "LINESTRING (0 0, 9 9)";
+    std::string b = "LINESTRING (4 2, 8 6)";
+    checkIntersectsDisjoint(a, b, false);
+    checkContainsWithin(a, b, false);
+}
+
+// testLinesClosedEmpty
+template<>
+template<>
+void object::test<17> ()
+{
+    std::string a = "MULTILINESTRING ((0 0, 0 1), (0 1, 1 1, 1 0, 0 0))";
+    std::string b = "LINESTRING EMPTY";
+    checkRelate(a, b, "FF1FFFFF2");
+    checkIntersectsDisjoint(a, b, false);
+    checkContainsWithin(a, b, false);
+}
+
+// testLinesRingTouchAtNode
+template<>
+template<>
+void object::test<18> ()
+{
+    std::string a = "LINESTRING (5 5, 1 8, 1 1, 5 5)";
+    std::string b = "LINESTRING (5 5, 9 5)";
+    checkRelate(a, b, "F01FFF102");
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, false);
+    checkTouches(a, b, true);
+}
+
+// testLinesTouchAtBdy
+template<>
+template<>
+void object::test<19> ()
+{
+    std::string a = "LINESTRING (5 5, 1 8)";
+    std::string b = "LINESTRING (5 5, 9 5)";
+    checkRelate(a, b, "FF1F00102");
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, false);
+    checkTouches(a, b, true);
+}
+
+// testLinesOverlapWithDisjointLine
+template<>
+template<>
+void object::test<20> ()
+{
+    std::string a = "LINESTRING (1 1, 9 9)";
+    std::string b = "MULTILINESTRING ((2 2, 8 8), (6 2, 8 4))";
+    checkRelate(a, b, "101FF0102");
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, false);
+    checkOverlaps(a, b, true);
+}
+
+// testLinesDisjointOverlappingEnvelopes
+template<>
+template<>
+void object::test<21> ()
+{
+    std::string a = "LINESTRING (60 0, 20 80, 100 80, 80 120, 40 140)";
+    std::string b = "LINESTRING (60 40, 140 40, 140 160, 0 160)";
+    checkRelate(a, b, "FF1FF0102");
+    checkIntersectsDisjoint(a, b, false);
+    checkContainsWithin(a, b, false);
+    checkTouches(a, b, false);
+}
+
+/**
+* Case from https://github.com/locationtech/jts/issues/270
+* Strictly, the lines cross, since their interiors intersect
+* according to the Orientation predicate.
+* However, the computation of the intersection point is
+* non-robust, and reports it as being equal to the endpoint
+* POINT (-10 0.0000000000000012)
+* For consistency the relate algorithm uses the intersection node topology.
+*/
+// testLinesCross_JTS270
+template<>
+template<>
+void object::test<22> ()
+{
+    std::string a = "LINESTRING (0 0, -10 0.0000000000000012)";
+    std::string b = "LINESTRING (-9.999143275740073 -0.1308959557133398, -10 0.0000000000001054)";
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, false);
+    checkCoversCoveredBy(a, b, false);
+    checkCrosses(a, b, false);
+    checkOverlaps(a, b, false);
+    checkTouches(a, b, true);
+}
+
+// testLinesContained_JTS396
+template<>
+template<>
+void object::test<23> ()
+{
+    std::string a = "LINESTRING (1 0, 0 2, 0 0, 2 2)";
+    std::string b = "LINESTRING (0 0, 2 2)";
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, true);
+    checkCoversCoveredBy(a, b, true);
+    checkCrosses(a, b, false);
+    checkOverlaps(a, b, false);
+    checkTouches(a, b, false);
+}
+
+
+/**
+* This case shows that lines must be self-noded,
+* so that node topology is constructed correctly
+* (at least for some predicates).
+*/
+// testLinesContainedWithSelfIntersection
+template<>
+template<>
+void object::test<24> ()
+{
+    std::string a = "LINESTRING (2 0, 0 2, 0 0, 2 2)";
+    std::string b = "LINESTRING (0 0, 2 2)";
+  //checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, true);
+    checkCoversCoveredBy(a, b, true);
+    checkCrosses(a, b, false);
+    checkOverlaps(a, b, false);
+    checkTouches(a, b, false);
+}
+
+// testLineContainedInRing
+template<>
+template<>
+void object::test<25> ()
+{
+    std::string a = "LINESTRING(60 60, 100 100, 140 60)";
+    std::string b = "LINESTRING(100 100, 180 20, 20 20, 100 100)";
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(b, a, true);
+    checkCoversCoveredBy(b, a, true);
+    checkCrosses(a, b, false);
+    checkOverlaps(a, b, false);
+    checkTouches(a, b, false);
+}
+
+// see https://github.com/libgeos/geos/issues/933
+// testLineLineProperIntersection
+template<>
+template<>
+void object::test<26> ()
+{
+    std::string a = "MULTILINESTRING ((0 0, 1 1), (0.5 0.5, 1 0.1, -1 0.1))";
+    std::string b = "LINESTRING (0 0, 1 1)";
+  //checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, true);
+    checkCoversCoveredBy(a, b, true);
+    checkCrosses(a, b, false);
+    checkOverlaps(a, b, false);
+    checkTouches(a, b, false);
+}
+
+// testLineSelfIntersectionCollinear
+template<>
+template<>
+void object::test<27> ()
+{
+    std::string a = "LINESTRING (9 6, 1 6, 1 0, 5 6, 9 6)";
+    std::string b = "LINESTRING (9 9, 3 1)";
+    checkRelate(a, b, "0F1FFF102");
+}
+
+//======= A/P  =============
+
+// testPolygonPointInside
+template<>
+template<>
+void object::test<28> ()
+{
+    std::string a = "POLYGON ((0 10, 10 10, 10 0, 0 0, 0 10))";
+    std::string b = "POINT (1 1)";
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, true);
+}
+
+// testPolygonPointOutside
+template<>
+template<>
+void object::test<29> ()
+{
+    std::string a = "POLYGON ((10 0, 0 0, 0 10, 10 0))";
+    std::string b = "POINT (8 8)";
+    checkIntersectsDisjoint(a, b, false);
+    checkContainsWithin(a, b, false);
+}
+
+// testPolygonPointInBoundary
+template<>
+template<>
+void object::test<30> ()
+{
+    std::string a = "POLYGON ((10 0, 0 0, 0 10, 10 0))";
+    std::string b = "POINT (1 0)";
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, false);
+    checkCoversCoveredBy(a, b, true);
+}
+
+// testAreaPointInExterior
+template<>
+template<>
+void object::test<31> ()
+{
+    std::string a = "POLYGON ((1 5, 5 5, 5 1, 1 1, 1 5))";
+    std::string b = "POINT (7 7)";
+    checkRelate(a, b, "FF2FF10F2");
+    checkIntersectsDisjoint(a, b, false);
+    checkContainsWithin(a, b, false);
+    checkCoversCoveredBy(a, b, false);
+    checkTouches(a, b, false);
+    checkOverlaps(a, b, false);
+}
+
+//======= A/L  =============
+
+
+// testAreaLineContainedAtLineVertex
+template<>
+template<>
+void object::test<32> ()
+{
+    std::string a = "POLYGON ((1 5, 5 5, 5 1, 1 1, 1 5))";
+    std::string b = "LINESTRING (2 3, 3 5, 4 3)";
+    checkIntersectsDisjoint(a, b, true);
+  //checkContainsWithin(a, b, true);
+  //checkCoversCoveredBy(a, b, true);
+    checkTouches(a, b, false);
+    checkOverlaps(a, b, false);
+}
+
+// testAreaLineTouchAtLineVertex
+template<>
+template<>
+void object::test<33> ()
+{
+    std::string a = "POLYGON ((1 5, 5 5, 5 1, 1 1, 1 5))";
+    std::string b = "LINESTRING (1 8, 3 5, 5 8)";
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, false);
+    checkCoversCoveredBy(a, b, false);
+    checkTouches(a, b, true);
+    checkOverlaps(a, b, false);
+}
+
+// testPolygonLineInside
+template<>
+template<>
+void object::test<34> ()
+{
+    std::string a = "POLYGON ((0 10, 10 10, 10 0, 0 0, 0 10))";
+    std::string b = "LINESTRING (1 8, 3 5, 5 8)";
+    checkRelate(a, b, "102FF1FF2");
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, true);
+}
+
+// testPolygonLineOutside
+template<>
+template<>
+void object::test<35> ()
+{
+    std::string a = "POLYGON ((10 0, 0 0, 0 10, 10 0))";
+    std::string b = "LINESTRING (4 8, 9 3)";
+    checkIntersectsDisjoint(a, b, false);
+    checkContainsWithin(a, b, false);
+}
+
+// testPolygonLineInBoundary
+template<>
+template<>
+void object::test<36> ()
+{
+    std::string a = "POLYGON ((10 0, 0 0, 0 10, 10 0))";
+    std::string b = "LINESTRING (1 0, 9 0)";
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, false);
+    checkCoversCoveredBy(a, b, true);
+    checkTouches(a, b, true);
+    checkOverlaps(a, b, false);
+}
+
+// testPolygonLineCrossingContained
+template<>
+template<>
+void object::test<37> ()
+{
+    std::string a = "MULTIPOLYGON (((20 80, 180 80, 100 0, 20 80)), ((20 160, 180 160, 100 80, 20 160)))";
+    std::string b = "LINESTRING (100 140, 100 40)";
+    checkRelate(a, b, "1020F1FF2");
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, true);
+    checkCoversCoveredBy(a, b, true);
+    checkTouches(a, b, false);
+    checkOverlaps(a, b, false);
+}
+
+// testValidateRelateLA_220
+template<>
+template<>
+void object::test<38> ()
+{
+    std::string a = "LINESTRING (90 210, 210 90)";
+    std::string b = "POLYGON ((150 150, 410 150, 280 20, 20 20, 150 150))";
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, false);
+    checkCoversCoveredBy(a, b, false);
+    checkTouches(a, b, false);
+    checkOverlaps(a, b, false);
+}
+
+  /**
+   * See RelateLA.xml (line 585)
+   */
+// testLineCrossingPolygonAtShellHolePoint
+template<>
+template<>
+void object::test<39> ()
+{
+    std::string a = "LINESTRING (60 160, 150 70)";
+    std::string b = "POLYGON ((190 190, 360 20, 20 20, 190 190), (110 110, 250 100, 140 30, 110 110))";
+    checkRelate(a, b, "F01FF0212");
+    checkTouches(a, b, true);
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, false);
+    checkCoversCoveredBy(a, b, false);
+    checkTouches(a, b, true);
+    checkOverlaps(a, b, false);
+}
+
+// testLineCrossingPolygonAtNonVertex
+template<>
+template<>
+void object::test<40> ()
+{
+    std::string a = "LINESTRING (20 60, 150 60)";
+    std::string b = "POLYGON ((150 150, 410 150, 280 20, 20 20, 150 150))";
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, false);
+    checkCoversCoveredBy(a, b, false);
+    checkTouches(a, b, false);
+    checkOverlaps(a, b, false);
+}
+
+// testPolygonLinesContainedCollinearEdge
+template<>
+template<>
+void object::test<41> ()
+{
+    std::string a = "POLYGON ((110 110, 200 20, 20 20, 110 110))";
+    std::string b = "MULTILINESTRING ((110 110, 60 40, 70 20, 150 20, 170 40), (180 30, 40 30, 110 80))";
+    checkRelate(a, b, "102101FF2");
+}
+
+//======= A/A  =============
+
+
+// testPolygonsEdgeAdjacent
+template<>
+template<>
+void object::test<42> ()
+{
+    std::string a = "POLYGON ((1 3, 3 3, 3 1, 1 1, 1 3))";
+    std::string b = "POLYGON ((5 3, 5 1, 3 1, 3 3, 5 3))";
+  //checkIntersectsDisjoint(a, b, true);
+    checkOverlaps(a, b, false);
+    checkTouches(a, b, true);
+    checkOverlaps(a, b, false);
+}
+
+// testPolygonsEdgeAdjacent2
+template<>
+template<>
+void object::test<43> ()
+{
+    std::string a = "POLYGON ((1 3, 4 3, 3 0, 1 1, 1 3))";
+    std::string b = "POLYGON ((5 3, 5 1, 3 0, 4 3, 5 3))";
+  //checkIntersectsDisjoint(a, b, true);
+    checkOverlaps(a, b, false);
+    checkTouches(a, b, true);
+    checkOverlaps(a, b, false);
+}
+
+// testPolygonsNested
+template<>
+template<>
+void object::test<44> ()
+{
+    std::string a = "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))";
+    std::string b = "POLYGON ((2 8, 8 8, 8 2, 2 2, 2 8))";
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, true);
+    checkCoversCoveredBy(a, b, true);
+    checkOverlaps(a, b, false);
+    checkTouches(a, b, false);
+}
+
+// testPolygonsOverlapProper
+template<>
+template<>
+void object::test<45> ()
+{
+    std::string a = "POLYGON ((1 1, 1 7, 7 7, 7 1, 1 1))";
+    std::string b = "POLYGON ((2 8, 8 8, 8 2, 2 2, 2 8))";
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, false);
+    checkCoversCoveredBy(a, b, false);
+    checkOverlaps(a, b, true);
+    checkTouches(a, b, false);
+}
+
+// testPolygonsOverlapAtNodes
+template<>
+template<>
+void object::test<46> ()
+{
+    std::string a = "POLYGON ((1 5, 5 5, 5 1, 1 1, 1 5))";
+    std::string b = "POLYGON ((7 3, 5 1, 3 3, 5 5, 7 3))";
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, false);
+    checkCoversCoveredBy(a, b, false);
+    checkOverlaps(a, b, true);
+    checkTouches(a, b, false);
+}
+
+// testPolygonsContainedAtNodes
+template<>
+template<>
+void object::test<47> ()
+{
+    std::string a = "POLYGON ((1 5, 5 5, 6 2, 1 1, 1 5))";
+    std::string b = "POLYGON ((1 1, 5 5, 6 2, 1 1))";
+  //checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, true);
+    checkCoversCoveredBy(a, b, true);
+    checkOverlaps(a, b, false);
+    checkTouches(a, b, false);
+}
+
+// testPolygonsNestedWithHole
+template<>
+template<>
+void object::test<48> ()
+{
+    std::string a = "POLYGON ((40 60, 420 60, 420 320, 40 320, 40 60), (200 140, 160 220, 260 200, 200 140))";
+    std::string b = "POLYGON ((80 100, 360 100, 360 280, 80 280, 80 100))";
+  //checkIntersectsDisjoint(true, a, b);
+    checkContainsWithin(a, b, false);
+    checkContainsWithin(b, a, false);
+  //checkCoversCoveredBy(false, a, b);
+  //checkOverlaps(true, a, b);
+    checkPredicate(*RelatePredicate::contains(), a, b, false);
+  //checkTouches(false, a, b);
+}
+
+// testPolygonsOverlappingWithBoundaryInside
+template<>
+template<>
+void object::test<49> ()
+{
+    std::string a = "POLYGON ((100 60, 140 100, 100 140, 60 100, 100 60))";
+    std::string b = "MULTIPOLYGON (((80 40, 120 40, 120 80, 80 80, 80 40)), ((120 80, 160 80, 160 120, 120 120, 120 80)), ((80 120, 120 120, 120 160, 80 160, 80 120)), ((40 80, 80 80, 80 120, 40 120, 40 80)))";
+    checkRelate(a, b, "21210F212");
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, false);
+    checkContainsWithin(b, a, false);
+    checkCoversCoveredBy(a, b, false);
+    checkOverlaps(a, b, true);
+    checkTouches(a, b, false);
+}
+
+// testPolygonsOverlapVeryNarrow
+template<>
+template<>
+void object::test<50> ()
+{
+    std::string a = "POLYGON ((120 100, 120 200, 200 200, 200 100, 120 100))";
+    std::string b = "POLYGON ((100 100, 100000 110, 100000 100, 100 100))";
+    checkRelate(a, b, "212111212");
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, false);
+    checkContainsWithin(b, a, false);
+  //checkCoversCoveredBy(false, a, b);
+  //checkOverlaps(true, a, b);
+  //checkTouches(false, a, b);
+}
+
+// testValidateRelateAA_86
+template<>
+template<>
+void object::test<51> ()
+{
+    std::string a = "POLYGON ((170 120, 300 120, 250 70, 120 70, 170 120))";
+    std::string b = "POLYGON ((150 150, 410 150, 280 20, 20 20, 150 150), (170 120, 330 120, 260 50, 100 50, 170 120))";
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, false);
+    checkCoversCoveredBy(a, b, false);
+    checkOverlaps(a, b, false);
+    checkPredicate(*RelatePredicate::within(), a, b, false);
+    checkTouches(a, b, true);
+}
+
+// testValidateRelateAA_97
+template<>
+template<>
+void object::test<52> ()
+{
+    std::string a = "POLYGON ((330 150, 200 110, 150 150, 280 190, 330 150))";
+    std::string b = "MULTIPOLYGON (((140 110, 260 110, 170 20, 50 20, 140 110)), ((300 270, 420 270, 340 190, 220 190, 300 270)))";
+    checkIntersectsDisjoint(a, b, true);
+    checkContainsWithin(a, b, false);
+    checkCoversCoveredBy(a, b, false);
+    checkOverlaps(a, b, false);
+    checkPredicate(*RelatePredicate::within(), a, b, false);
+    checkTouches(a, b, true);
+}
+
+// testAdjacentPolygons
+template<>
+template<>
+void object::test<53> ()
+{
+    std::string a = "POLYGON ((1 9, 6 9, 6 1, 1 1, 1 9))";
+    std::string b = "POLYGON ((9 9, 9 4, 6 4, 6 9, 9 9))";
+    checkRelateMatches(a, b, IntersectionMatrixPattern::ADJACENT, true);
+}
+
+// testAdjacentPolygonsTouchingAtPoint
+template<>
+template<>
+void object::test<54> ()
+{
+    std::string a = "POLYGON ((1 9, 6 9, 6 1, 1 1, 1 9))";
+    std::string b = "POLYGON ((9 9, 9 4, 6 4, 7 9, 9 9))";
+    checkRelateMatches(a, b, IntersectionMatrixPattern::ADJACENT, false);
+}
+
+// testAdjacentPolygonsOverlappping
+template<>
+template<>
+void object::test<55> ()
+{
+    std::string a = "POLYGON ((1 9, 6 9, 6 1, 1 1, 1 9))";
+    std::string b = "POLYGON ((9 9, 9 4, 6 4, 5 9, 9 9))";
+    checkRelateMatches(a, b, IntersectionMatrixPattern::ADJACENT, false);
+}
+
+// testContainsProperlyPolygonContained
+template<>
+template<>
+void object::test<56> ()
+{
+    std::string a = "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))";
+    std::string b = "POLYGON ((2 8, 5 8, 5 5, 2 5, 2 8))";
+    checkRelateMatches(a, b, IntersectionMatrixPattern::CONTAINS_PROPERLY, true);
+}
+
+// testContainsProperlyPolygonTouching
+template<>
+template<>
+void object::test<57> ()
+{
+    std::string a = "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))";
+    std::string b = "POLYGON ((9 1, 5 1, 5 5, 9 5, 9 1))";
+    checkRelateMatches(a, b, IntersectionMatrixPattern::CONTAINS_PROPERLY, false);
+}
+
+// testContainsProperlyPolygonsOverlapping
+template<>
+template<>
+void object::test<58> ()
+{
+    std::string a = "GEOMETRYCOLLECTION (POLYGON ((1 9, 6 9, 6 4, 1 4, 1 9)), POLYGON ((2 4, 6 7, 9 1, 2 4)))";
+    std::string b = "POLYGON ((5 5, 6 5, 6 4, 5 4, 5 5))";
+    checkRelateMatches(a, b, IntersectionMatrixPattern::CONTAINS_PROPERLY, true);
+}
+
+//================  Repeated Points  ==============
+
+// testRepeatedPointLL
+template<>
+template<>
+void object::test<59> ()
+{
+    std::string a = "LINESTRING(0 0, 5 5, 5 5, 5 5, 9 9)";
+    std::string b = "LINESTRING(0 9, 5 5, 5 5, 5 5, 9 0)";
+    checkRelate(a, b, "0F1FF0102");
+    checkIntersectsDisjoint(a, b, true);
+}
+
+// testRepeatedPointAA
+template<>
+template<>
+void object::test<60> ()
+{
+    std::string a = "POLYGON ((1 9, 9 7, 9 1, 1 3, 1 9))";
+    std::string b = "POLYGON ((1 3, 1 3, 1 3, 3 7, 9 7, 9 7, 1 3))";
+    checkRelate(a, b, "212F01FF2");
+}
+
+
+
+
+
+
+} // namespace tut

commit 5b8e2679e2b0db5d458007abe392ee0a8d8fd8d2
Merge: 8ee1ae84d 4907d593a
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date:   Thu Jul 25 11:41:49 2024 -0700

    Merge branch 'main' of github.com:libgeos/geos into main-relate-ng


-----------------------------------------------------------------------

Summary of changes:
 include/geos/geom/GeometryCollection.h             |   7 -
 .../util/{GeometryExtracter.h => GeometryLister.h} |  45 +-
 .../noding/MCIndexSegmentSetMutualIntersector.h    |  21 +-
 .../operation/relateng/EdgeSegmentIntersector.h    |   4 +-
 .../geos/operation/relateng/EdgeSetIntersector.h   |   8 +-
 .../operation/relateng/IntersectionMatrixPattern.h |   6 +-
 include/geos/operation/relateng/NodeSections.h     |   4 +-
 include/geos/operation/relateng/RelateGeometry.h   |  40 +-
 include/geos/operation/relateng/RelateNG.h         | 261 ++++++
 .../geos/operation/relateng/RelateSegmentString.h  |  31 +-
 include/geos/operation/relateng/TopologyComputer.h |  36 +-
 .../geos/operation/relateng/TopologyPredicate.h    |   7 +
 src/geom/GeometryCollection.cpp                    |  17 -
 src/noding/MCIndexSegmentSetMutualIntersector.cpp  |  13 +-
 src/operation/relateng/EdgeSetIntersector.cpp      |   8 +-
 src/operation/relateng/RelateGeometry.cpp          |  96 ++-
 src/operation/relateng/RelateNG.cpp                | 523 ++++++++++++
 src/operation/relateng/RelateSegmentString.cpp     |  57 +-
 src/operation/relateng/TopologyComputer.cpp        |  79 +-
 tests/unit/geom/util/GeometryListerTest.cpp        |  62 ++
 tests/unit/operation/relateng/RelateNGTest.cpp     | 935 +++++++++++++++++++++
 21 files changed, 2041 insertions(+), 219 deletions(-)
 copy include/geos/geom/util/{GeometryExtracter.h => GeometryLister.h} (55%)
 create mode 100644 include/geos/operation/relateng/RelateNG.h
 create mode 100644 src/operation/relateng/RelateNG.cpp
 create mode 100644 tests/unit/geom/util/GeometryListerTest.cpp
 create mode 100644 tests/unit/operation/relateng/RelateNGTest.cpp


hooks/post-receive
-- 
GEOS


More information about the geos-commits mailing list