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

[git at osgeo.org]

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  4faaea0f46da942c0e3a0db9f6a5f96aa0cb8e17 (commit)
       via  cc2ad9af20c5a2ff050afae96518c19a79fd5ae3 (commit)
       via  23b7399ba7510ab7a2462c2c4e47867b32091f97 (commit)
       via  6201229a66ae3eef2b0c62cc39dbab48b2e3df7e (commit)
       via  42546119c35e65aad72dea1477eb4a057ead631e (commit)
       via  8e525c825caf314d61a6279a426a59ded115681e (commit)
       via  d3b22f1393eb31279647b1fe2aadc1e804a15fe2 (commit)
       via  3db005a65baf5fb53a3d2078c30dd7d365115b5e (commit)
      from  65803b8f139ffa4d4c36aa0811a3827d3fa3d280 (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 4faaea0f46da942c0e3a0db9f6a5f96aa0cb8e17
Merge: cc2ad9af2 6201229a6
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date:   Mon Jul 22 15:27:03 2024 -0700

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


commit cc2ad9af20c5a2ff050afae96518c19a79fd5ae3
Merge: 23b7399ba 65803b8f1
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date:   Mon Jul 22 15:26:50 2024 -0700

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


commit 23b7399ba7510ab7a2462c2c4e47867b32091f97
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date:   Mon Jul 22 15:24:34 2024 -0700

    Port RelateSegmentString, RelateGeometry

diff --git a/include/geos/geom/GeometryCollection.h b/include/geos/geom/GeometryCollection.h
index ae9a441aa..27de4c7cf 100644
--- a/include/geos/geom/GeometryCollection.h
+++ b/include/geos/geom/GeometryCollection.h
@@ -196,6 +196,13 @@ 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);
@@ -248,6 +255,7 @@ protected:
 
     bool hasCurvedComponents() const override;
 
+
 };
 
 } // namespace geos::geom
diff --git a/include/geos/geom/IntersectionMatrix.h b/include/geos/geom/IntersectionMatrix.h
index 76823abbd..b2078ff5a 100644
--- a/include/geos/geom/IntersectionMatrix.h
+++ b/include/geos/geom/IntersectionMatrix.h
@@ -226,7 +226,7 @@ public:
      *
      * @return the dimension value at the given matrix position.
      */
-    int get(geom::Location row, geom::Location column) const {
+    int get(Location row, Location column) const {
         return matrix[static_cast<size_t>(row)][static_cast<size_t>(column)];
     }
 
@@ -363,6 +363,7 @@ public:
      */
     std::string toString() const;
 
+
 private:
 
     static const int firstDim; // = 3;
diff --git a/include/geos/noding/SegmentString.h b/include/geos/noding/SegmentString.h
index 3d71f009a..af7cee2de 100644
--- a/include/geos/noding/SegmentString.h
+++ b/include/geos/noding/SegmentString.h
@@ -131,6 +131,42 @@ public:
         return ss.getSegmentOctant(index);
     }
 
+    /**
+     * Gets the next vertex in a ring from a vertex index.
+     *
+     * @param index the vertex index
+     * @return the next vertex in the ring
+     *
+     * @see isClosed
+     */
+    const geom::CoordinateXY& nextInRing(std::size_t index) const
+    {
+        std::size_t nextIndex = index + 1;
+        if (nextIndex > size() - 1) {
+            nextIndex = 1;
+        }
+        return getCoordinate(nextIndex);
+    }
+
+    /**
+     * Gets the previous vertex in a ring from a vertex index.
+     *
+     * @param index the vertex index
+     * @return the previous vertex in the ring
+     *
+     * @see isClosed
+     */
+    const geom::CoordinateXY& prevInRing(std::size_t index) const
+    {
+        std::size_t prevIndex;
+        if (index == 0)
+            prevIndex = size() - 2;
+        else
+            prevIndex = index - 1;
+        return getCoordinate( prevIndex );
+    }
+
+
     bool isClosed() const {
         return seq->front<geom::CoordinateXY>().equals(seq->back<geom::CoordinateXY>());
     }
diff --git a/include/geos/operation/relateng/AdjacentEdgeLocator.h b/include/geos/operation/relateng/AdjacentEdgeLocator.h
index fe72d5df9..d6e576985 100644
--- a/include/geos/operation/relateng/AdjacentEdgeLocator.h
+++ b/include/geos/operation/relateng/AdjacentEdgeLocator.h
@@ -15,6 +15,7 @@
 
 #pragma once
 
+#include <geos/geom/CoordinateSequence.h>
 #include <geos/geom/Location.h>
 #include <geos/export.h>
 
@@ -28,7 +29,6 @@ namespace relateng {
 }
 }
 namespace geom {
-    class CoordinateSequence;
     class CoordinateXY;
     class Geometry;
     class LinearRing;
diff --git a/include/geos/operation/relateng/BasicPredicate.h b/include/geos/operation/relateng/BasicPredicate.h
new file mode 100644
index 000000000..2ed3f0539
--- /dev/null
+++ b/include/geos/operation/relateng/BasicPredicate.h
@@ -0,0 +1,105 @@
+/**********************************************************************
+ *
+ * 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/geom/Location.h>
+#include <geos/operation/relateng/TopologyPredicate.h>
+#include <geos/export.h>
+
+#include <memory>
+#include <vector>
+
+// Forward declarations
+namespace geos {
+namespace geom {
+    class Envelope;
+}
+}
+
+
+using geos::geom::Envelope;
+using geos::geom::Location;
+
+
+namespace geos {      // geos.
+namespace operation { // geos.operation.
+namespace relateng {  // geos.operation.relateng
+
+
+class GEOS_DLL BasicPredicate : public TopologyPredicate {
+
+private:
+
+    static constexpr int UNKNOWN = -1;
+    static constexpr int FALSE = 0;
+    static constexpr int TRUE = 1;
+
+    int m_value = UNKNOWN;
+
+    static bool isKnown(int val);
+
+    static bool toBoolean(int val);
+
+    static int toValue(bool val);
+
+
+protected:
+
+    /**
+    * Updates the predicate value to the given state
+    * if it is currently unknown.
+    *
+    * @param val the predicate value to update
+    */
+    void setValue(bool val);
+
+    void setValue(int val);
+
+    void setValueIf(bool val, bool cond);
+
+    void require(bool cond);
+
+    using TopologyPredicate::requireCovers;
+    void requireCovers(const Envelope& a, const Envelope& b);
+
+
+public:
+
+    /**
+    * Tests if two geometries intersect
+    * based on an interaction at given locations.
+    *
+    * @param locA the location on geometry A
+    * @param locB the location on geometry B
+    * @return true if the geometries intersect
+    */
+    static bool isIntersection(Location locA, Location locB);
+
+    std::string name() const override = 0;
+
+    void finish() override = 0;
+
+    bool isKnown() const override;
+
+    bool value() const override;
+
+
+};
+
+} // namespace geos.operation.relateng
+} // namespace geos.operation
+} // namespace geos
+
diff --git a/include/geos/operation/relateng/IMPatternMatcher.h b/include/geos/operation/relateng/IMPatternMatcher.h
new file mode 100644
index 000000000..be8486c57
--- /dev/null
+++ b/include/geos/operation/relateng/IMPatternMatcher.h
@@ -0,0 +1,87 @@
+/**********************************************************************
+ *
+ * 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/geom/Location.h>
+#include <geos/geom/Dimension.h>
+#include <geos/geom/IntersectionMatrix.h>
+#include <geos/operation/relateng/IMPredicate.h>
+#include <geos/export.h>
+
+#include <memory>
+#include <vector>
+
+// Forward declarations
+namespace geos {
+namespace geom {
+    class Envelope;
+}
+}
+
+
+using geos::geom::Envelope;
+using geos::geom::Location;
+using geos::geom::Dimension;
+using geos::geom::IntersectionMatrix;
+
+
+namespace geos {      // geos.
+namespace operation { // geos.operation.
+namespace relateng {  // geos.operation.relateng
+
+
+class GEOS_DLL IMPatternMatcher : public IMPredicate {
+
+
+private:
+
+    std::string imPattern;
+    IntersectionMatrix patternMatrix;
+
+    static bool requireInteraction(const IntersectionMatrix& im);
+
+    static bool isInteraction(int imDim);
+
+
+public:
+
+    IMPatternMatcher(std::string p_imPattern)
+        : imPattern(p_imPattern)
+        , patternMatrix(p_imPattern)
+        {};
+
+    std::string name() const override;
+
+    using IMPredicate::init;
+    void init(const Envelope& envA, const Envelope& envB) override;
+
+    bool requireInteraction() const override;
+
+    bool isDetermined() const override;
+
+    bool valueIM() override;
+
+    std::string toString() const;
+
+    friend std::ostream& operator<<(std::ostream& os, const IMPatternMatcher& imp);
+
+};
+
+
+} // namespace geos.operation.relateng
+} // namespace geos.operation
+} // namespace geos
+
diff --git a/include/geos/operation/relateng/IMPredicate.h b/include/geos/operation/relateng/IMPredicate.h
new file mode 100644
index 000000000..6cf9f9b10
--- /dev/null
+++ b/include/geos/operation/relateng/IMPredicate.h
@@ -0,0 +1,131 @@
+/**********************************************************************
+ *
+ * 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/geom/Location.h>
+#include <geos/geom/Dimension.h>
+#include <geos/geom/IntersectionMatrix.h>
+#include <geos/operation/relateng/BasicPredicate.h>
+#include <geos/export.h>
+
+#include <memory>
+#include <vector>
+
+// Forward declarations
+namespace geos {
+namespace geom {
+    class Envelope;
+}
+}
+
+
+using geos::geom::Envelope;
+using geos::geom::Location;
+using geos::geom::Dimension;
+using geos::geom::IntersectionMatrix;
+
+
+namespace geos {      // geos.
+namespace operation { // geos.operation.
+namespace relateng {  // geos.operation.relateng
+
+
+class GEOS_DLL IMPredicate : public BasicPredicate {
+
+private:
+
+
+
+protected:
+
+    static constexpr int DIM_UNKNOWN = Dimension::DONTCARE;
+
+    int dimA;
+    int dimB;
+    IntersectionMatrix intMatrix;
+
+    /**
+     * Gets the value of the predicate according to the current
+     * intersection matrix state.
+     *
+     * @return the current predicate value
+     */
+    virtual bool valueIM() = 0;
+
+    /**
+     * Tests whether predicate evaluation can be short-circuited
+     * due to the current state of the matrix providing
+     * enough information to determine the predicate value.
+     *
+     * If this value is true then valueIM()
+     * must provide the correct result of the predicate.
+     *
+     * @return true if the predicate value is determined
+     */
+    virtual bool isDetermined() const = 0;
+
+    /**
+     * Tests whether the exterior of the specified input geometry
+     * is intersected by any part of the other input.
+     *
+     * @param isA the input geometry
+     * @return true if the input geometry exterior is intersected
+     */
+    bool intersectsExteriorOf(bool isA) const;
+
+    bool isIntersects(Location locA, Location locB) const;
+
+
+public:
+
+    IMPredicate()
+    {
+        // intMatrix = new IntersectionMatrix();
+        //-- E/E is always dim = 2
+        intMatrix.set(Location::EXTERIOR, Location::EXTERIOR, Dimension::A);
+    }
+
+    static bool isDimsCompatibleWithCovers(int dim0, int dim1);
+
+    void init(int dA, int dB) override;
+
+    void updateDimension(Location locA, Location locB, int dimension) override;
+
+    bool isDimChanged(Location locA, Location locB, int dimension) const;
+
+    using TopologyPredicate::isKnown;
+    bool isKnown(Location locA, Location locB) const;
+
+    bool isDimension(Location locA, Location locB, int dimension) const;
+
+    int getDimension(Location locA, Location locB) const;
+
+    /**
+     * Sets the final value based on the state of the IM.
+     */
+    void finish() override;
+
+    std::string toString() const;
+
+    friend std::ostream& operator<<(std::ostream& os, const IMPredicate& imp);
+
+
+};
+
+} // namespace geos.operation.relateng
+} // namespace geos.operation
+} // namespace geos
+
diff --git a/include/geos/operation/relateng/IntersectionMatrixPattern.h b/include/geos/operation/relateng/IntersectionMatrixPattern.h
new file mode 100644
index 000000000..a7ef3462c
--- /dev/null
+++ b/include/geos/operation/relateng/IntersectionMatrixPattern.h
@@ -0,0 +1,65 @@
+/**********************************************************************
+ *
+ * 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 <string>
+
+
+namespace geos {      // geos.
+namespace operation { // geos.operation
+namespace relateng {  // geos.operation.relateng
+
+
+class GEOS_DLL IntersectionMatrixPattern {
+
+private:
+
+    /**
+     * Cannot be instantiated.
+     */
+    IntersectionMatrixPattern() {};
+
+
+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****";
+
+    /**
+     * 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*";
+
+    /**
+     * 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********";
+
+
+};
+
+} // namespace geos.operation.relateng
+} // namespace geos.operation
+} // namespace geos
+
diff --git a/include/geos/operation/relateng/RelateGeometry.h b/include/geos/operation/relateng/RelateGeometry.h
index 3969f414e..59ebddcbb 100644
--- a/include/geos/operation/relateng/RelateGeometry.h
+++ b/include/geos/operation/relateng/RelateGeometry.h
@@ -15,22 +15,39 @@
 
 #pragma once
 
-//#include <geos/operation/relateng/NodeSection.h>
+#include <geos/algorithm/BoundaryNodeRule.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/Dimension.h>
+#include <geos/geom/Location.h>
+#include <geos/operation/relateng/RelatePointLocator.h>
 #include <geos/export.h>
 
 
+#include <string>
+#include <sstream>
+
+
 // Forward declarations
 namespace geos {
+namespace geom {
+    class CoordinateSequence;
+    class Envelope;
+    class Geometry;
+    class LinearRing;
+    class LineString;
+    class MultiPolygon;
+    class Point;
+}
 namespace operation {
 namespace relateng {
-// class RelateNode;
+    class RelateSegmentString;
 }
 }
 }
 
 
-// using geos::geom::CoordinateXY;
-// using geos::geom::Geometry;
+using geos::algorithm::BoundaryNodeRule;
+using namespace geos::geom;
 
 
 namespace geos {      // geos.
@@ -44,14 +61,150 @@ private:
 
     // Members
 
+    const Geometry* geom;
+    bool m_isPrepared = false;
+    const Envelope* geomEnv;
+    const BoundaryNodeRule& boundaryNodeRule;
+    int geomDim = Dimension::False;
+    bool isLineZeroLen = false;
+    bool isGeomEmpty = false;
+
+    Coordinate::ConstXYSet uniquePoints;
+    std::unique_ptr<RelatePointLocator> locator;
+    int elementId = 0;
+    bool hasPoints = false;
+    bool hasLines = false;
+    bool hasAreas = false;
+
+
     // Methods
 
+    void analyzeDimensions();
+
+    /**
+    * Tests if all geometry linear elements are zero-length.
+    * For efficiency the test avoids computing actual length.
+    *
+    * @param geom
+    * @return
+    */
+    static bool isZeroLength(const Geometry* geom);
+
+    static bool isZeroLength(const LineString* line);
+
+    RelatePointLocator* getLocator();
+
+    Coordinate::ConstXYSet createUniquePoints();
+
+    void extractSegmentStringsFromAtomic(bool isA,
+        const Geometry* geom, const MultiPolygon* parentPolygonal,
+        const Envelope* env,
+        std::vector<std::unique_ptr<RelateSegmentString>>& segStrings);
+
+    void extractRingToSegmentString(bool isA,
+        const LinearRing* ring, int ringId, const Envelope* env,
+        const Geometry* parentPoly,
+        std::vector<std::unique_ptr<RelateSegmentString>>& segStrings);
+
+    void extractSegmentStrings(bool isA,
+        const Envelope* env, const Geometry* geom,
+        std::vector<std::unique_ptr<RelateSegmentString>>& segStrings);
+
 
 public:
 
     static constexpr bool GEOM_A = true;
     static constexpr bool GEOM_B = false;
 
+    RelateGeometry(const Geometry* input)
+        : RelateGeometry(input, false, BoundaryNodeRule::getBoundaryRuleMod2())
+        {};
+
+    RelateGeometry(const Geometry* input, const BoundaryNodeRule& bnRule)
+        : RelateGeometry(input, false, bnRule)
+        {};
+
+    RelateGeometry(const Geometry* input, bool p_isPrepared, const BoundaryNodeRule& bnRule);
+
+    static std::string name(bool isA);
+
+    const Geometry* getGeometry() const;
+
+    bool isPrepared() const;
+
+    const Envelope* getEnvelope() const;
+
+    int getDimension() const;
+
+    bool hasDimension(int dim) const;
+
+    /**
+    * Gets the actual non-empty dimension of the geometry.
+    * Zero-length LineStrings are treated as Points.
+    *
+    * @return the real (non-empty) dimension
+    */
+    int getDimensionReal() const;
+
+    bool hasEdges() const;
+
+    bool isNodeInArea(const CoordinateXY* nodePt, const Geometry* parentPolygonal);
+
+    Location locateLineEnd(const CoordinateXY* p) ;
+
+    /**
+     * Locates a vertex of a polygon.
+     * A vertex of a Polygon or MultiPolygon is on
+     * the {@link Location#BOUNDARY}.
+     * But a vertex of an overlapped polygon in a GeometryCollection
+     * may be in the {@link Location#INTERIOR}.
+     *
+     * @param pt the polygon vertex
+     * @return the location of the vertex
+     */
+    Location locateAreaVertex(const CoordinateXY* pt);
+
+    Location locateNode(const CoordinateXY* pt, const Geometry* parentPolygonal);
+
+    int locateWithDim(const CoordinateXY* pt);
+
+    bool isPointsOrPolygons() const;
+
+    /**
+     * Tests whether the geometry has polygonal topology.
+     * This is not the case if it is a GeometryCollection
+     * containing more than one polygon (since they may overlap
+     * or be adjacent).
+     * The significance is that polygonal topology allows more assumptions
+     * about the location of boundary vertices.
+     *
+     * @return true if the geometry has polygonal topology
+     */
+    bool isPolygonal() const;
+
+    bool isEmpty() const;
+
+    bool hasBoundary();
+
+    Coordinate::ConstXYSet& getUniquePoints();
+
+    std::vector<const Point*> getEffectivePoints();
+
+    /**
+     * Extract RelateSegmentStrings 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
+     */
+    std::vector<std::unique_ptr<RelateSegmentString>> extractSegmentStrings(bool isA, const Envelope* env);
+
+    std::string toString() const;
+
+    friend std::ostream& operator<<(std::ostream& os, const RelateGeometry& rg);
+
 };
 
 } // namespace geos.operation.relateng
diff --git a/include/geos/operation/relateng/RelateMatrixPredicate.h b/include/geos/operation/relateng/RelateMatrixPredicate.h
new file mode 100644
index 000000000..425884cd5
--- /dev/null
+++ b/include/geos/operation/relateng/RelateMatrixPredicate.h
@@ -0,0 +1,85 @@
+/**********************************************************************
+ *
+ * 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/geom/Location.h>
+#include <geos/geom/Dimension.h>
+#include <geos/geom/IntersectionMatrix.h>
+#include <geos/operation/relateng/IMPredicate.h>
+#include <geos/export.h>
+
+#include <memory>
+#include <vector>
+
+// Forward declarations
+namespace geos {
+namespace geom {
+    class Envelope;
+}
+}
+
+
+using geos::geom::Envelope;
+using geos::geom::Location;
+using geos::geom::Dimension;
+using geos::geom::IntersectionMatrix;
+
+
+namespace geos {      // geos.
+namespace operation { // geos.operation.
+namespace relateng {  // geos.operation.relateng
+
+
+class GEOS_DLL RelateMatrixPredicate : public IMPredicate {
+
+public:
+
+    RelateMatrixPredicate() {};
+
+    std::string name() const override {
+        return "relateMatrix";
+    };
+
+    bool requireInteraction() const override {
+        //-- ensure entire matrix is computed
+        return false;
+    };
+
+    bool isDetermined() const override {
+        //-- ensure entire matrix is computed
+        return false;
+    };
+
+    bool valueIM() override {
+        //-- indicates full matrix is being evaluated
+        return false;
+    };
+
+    /**
+    * Gets the current state of the IM matrix (which may only be partially complete).
+    *
+    * @return the IM matrix
+    */
+    IntersectionMatrix& getIM() {
+        return intMatrix;
+    }
+
+};
+
+} // namespace geos.operation.relateng
+} // namespace geos.operation
+} // namespace geos
+
diff --git a/include/geos/operation/relateng/RelatePredicate.h b/include/geos/operation/relateng/RelatePredicate.h
new file mode 100644
index 000000000..38a6cd676
--- /dev/null
+++ b/include/geos/operation/relateng/RelatePredicate.h
@@ -0,0 +1,639 @@
+/**********************************************************************
+ *
+ * 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/geom/Location.h>
+#include <geos/operation/relateng/BasicPredicate.h>
+#include <geos/operation/relateng/IMPatternMatcher.h>
+#include <geos/operation/relateng/IMPredicate.h>
+#include <geos/operation/relateng/RelateGeometry.h>
+#include <geos/geom/Envelope.h>
+#include <geos/export.h>
+
+
+using geos::geom::Envelope;
+using geos::geom::Location;
+
+
+namespace geos {      // geos.
+namespace operation { // geos.operation.
+namespace relateng {  // geos.operation.relateng
+
+
+class GEOS_DLL RelatePredicate {
+
+public:
+
+/************************************************************************
+ *
+ * Creates a predicate to determine whether two geometries intersect.
+ *
+ * The intersects predicate has the following equivalent definitions:
+ *
+ *  * The two geometries have at least one point in common
+ *  * The DE-9IM Intersection Matrix for the two geometries matches
+ *    at least one of the patterns
+ *
+ *    [T********]
+ *    [*T*******]
+ *    [***T*****]
+ *    [****T****]
+ *
+ *  disjoint() = false
+ *  (intersects is the inverse of disjoint)
+ *
+ * @return the predicate instance
+ *
+ * @see disjoint()
+ */
+class IntersectsPredicate : public BasicPredicate {
+
+public:
+
+    std::string name() const override {
+        return std::string("intersects");
+    }
+
+    bool requireSelfNoding() const override {
+        //-- self-noding is not required to check for a simple interaction
+        return false;
+    }
+
+    bool requireExteriorCheck(bool isSourceA) const override {
+        (void)isSourceA;
+        //-- intersects only requires testing interaction
+        return false;
+    }
+
+    void init(const Envelope& envA, const Envelope& envB) override {
+        require(envA.intersects(envB));
+    }
+
+    void updateDimension(Location locA, Location locB, int dimension) override {
+        (void)dimension;
+        setValueIf(true, isIntersection(locA, locB));
+    }
+
+    void finish() override {
+        //-- if no intersecting locations were found
+        setValue(false);
+    }
+
+};
+
+static std::unique_ptr<BasicPredicate> intersects();
+
+/************************************************************************
+ *
+ * Creates a predicate to determine whether two geometries are disjoint.
+ *
+ * The disjoint predicate has the following equivalent definitions:
+ *
+ *   * The two geometries have no point in common
+ *   * The DE-9IM Intersection Matrix for the two geometries matches
+ *        [FF*FF****]
+ *   * intersects() = false
+ *     (disjoint is the inverse of intersects)
+ *
+ * @return the predicate instance
+ *
+ * @see intersects()
+ */
+class DisjointPredicate : public BasicPredicate {
+
+    std::string name() const override {
+        return std::string("disjoint");
+    }
+
+    bool requireSelfNoding() const override {
+        //-- self-noding is not required to check for a simple interaction
+        return false;
+    }
+
+    bool requireExteriorCheck(bool isSourceA) const override {
+        (void)isSourceA;
+        //-- intersects only requires testing interaction
+        return false;
+    }
+
+    void init(const Envelope& envA, const Envelope& envB) override {
+        setValueIf(true, envA.disjoint(envB));
+    }
+
+    void updateDimension(Location locA, Location locB, int dimension) override {
+        (void)dimension;
+        setValueIf(false, isIntersection(locA, locB));
+    }
+
+    void finish() override {
+        //-- if no intersecting locations were found
+        setValue(true);
+    }
+};
+
+static std::unique_ptr<BasicPredicate> disjoint();
+
+/************************************************************************
+ * Creates a predicate to determine whether a geometry contains another geometry.
+ *
+ * The contains predicate has the following equivalent definitions:
+ *
+ *   * Every point of the other geometry is a point of this geometry,
+ *     and the interiors of the two geometries have at least one point in common.
+ *   * The DE-9IM Intersection Matrix for the two geometries matches
+ *     the pattern
+ *       [T*****FF*]
+ *   * within(B, A) = true
+ *     (contains is the converse of within)
+ *
+ * An implication of the definition is that "Geometries do not
+ * contain their boundary".  In other words, if a geometry A is a subset of
+ * the points in the boundary of a geometry B, B.contains(A) = false.
+ * (As a concrete example, take A to be a LineString which lies in the boundary of a Polygon B.)
+ * For a predicate with similar behavior but avoiding
+ * this subtle limitation, see covers().
+ *
+ * @return the predicate instance
+ *
+ * @see within()
+ */
+class ContainsPredicate : public IMPredicate {
+
+    std::string name() const override {
+        return std::string("contains");
+    }
+
+    bool requireCovers(bool isSourceA) override {
+        return isSourceA == RelateGeometry::GEOM_A;
+    }
+
+    bool requireExteriorCheck(bool isSourceA) const override {
+        //-- only need to check B against Exterior of A
+        return isSourceA == RelateGeometry::GEOM_B;
+    }
+
+    void init(int dA, int dB) override {
+        IMPredicate::init(dA, dB);
+        require(isDimsCompatibleWithCovers(dA, dB));
+    }
+
+    void init(const Envelope& envA, const Envelope& envB) override {
+        BasicPredicate::requireCovers(envA, envB);
+    }
+
+    bool isDetermined() const override {
+        return intersectsExteriorOf(RelateGeometry::GEOM_A);
+    }
+
+    bool valueIM() override {
+        return intMatrix.isContains();
+    }
+};
+
+static std::unique_ptr<IMPredicate> contains();
+
+
+
+/************************************************************************
+ * Creates a predicate to determine whether a geometry is within another geometry.
+ *
+ * The within predicate has the following equivalent definitions:
+ *
+ *   * Every point of this geometry is a point of the other geometry,
+ *     and the interiors of the two geometries have at least one point in common.
+ *   * The DE-9IM Intersection Matrix for the two geometries matches
+ *     [T*F**F***]
+ *   *  contains(B, A) = true
+ *      (within is the converse of  contains())
+ *
+ * An implication of the definition is that
+ * "The boundary of a Geometry is not within the Geometry".
+ * In other words, if a geometry A is a subset of
+ * the points in the boundary of a geometry B, within(B, A) = false
+ * (As a concrete example, take A to be a LineString which lies in the boundary of a Polygon B.)
+ * For a predicate with similar behavior but avoiding
+ * this subtle limitation, see coveredBy().
+ *
+ * @return the predicate instance
+ *
+ * @see #contains()
+ */
+class WithinPredicate : public IMPredicate {
+
+    std::string name() const override {
+        return std::string("within");
+    }
+
+    bool requireCovers(bool isSourceA) override {
+        return isSourceA == RelateGeometry::GEOM_B;
+    }
+
+    bool requireExteriorCheck(bool isSourceA) const override {
+        //-- only need to check B against Exterior of A
+        return isSourceA == RelateGeometry::GEOM_A;
+    }
+
+    void init(int dA, int dB) override {
+        IMPredicate::init(dA, dB);
+        require(isDimsCompatibleWithCovers(dB, dA));
+    }
+
+    void init(const Envelope& envA, const Envelope& envB) override {
+        BasicPredicate::requireCovers(envB, envA);
+    }
+
+    bool isDetermined() const override {
+        return intersectsExteriorOf(RelateGeometry::GEOM_B);
+    }
+
+    bool valueIM() override {
+        return intMatrix.isWithin();
+    }
+};
+
+static std::unique_ptr<IMPredicate> within();
+
+
+
+/************************************************************************
+ * Creates a predicate to determine whether a geometry covers another geometry.
+ *
+ * The covers predicate has the following equivalent definitions:
+ *
+ * Every point of the other geometry is a point of this geometry.
+ * The DE-9IM Intersection Matrix for the two geometries matches
+ * at least one of the following patterns:
+ *
+ *  * [T*****FF*]
+ *  * [*T****FF*]
+ *  * [***T**FF*]
+ *  * [****T*FF*]
+ *
+ * coveredBy(b, a) = true
+ * (covers is the converse of coveredBy())
+ *
+ * If either geometry is empty, the value of this predicate is false.
+ *
+ * This predicate is similar to contains(),
+ * but is more inclusive (i.e. returns true for more cases).
+ * In particular, unlike contains it does not distinguish between
+ * points in the boundary and in the interior of geometries.
+ * For most cases, covers should be used in preference to contains.
+ * As an added benefit, covers is more amenable to optimization,
+ * and hence should be more performant.
+ *
+ * @return the predicate instance
+ *
+ * @see #coveredBy()
+ */
+class CoversPredicate : public IMPredicate {
+
+    std::string name() const override {
+        return std::string("covers");
+    }
+
+    bool requireCovers(bool isSourceA) override {
+        return isSourceA == RelateGeometry::GEOM_A;
+    }
+
+    bool requireExteriorCheck(bool isSourceA) const override {
+        //-- only need to check B against Exterior of A
+        return isSourceA == RelateGeometry::GEOM_B;
+    }
+
+    void init(int dA, int dB) override {
+        IMPredicate::init(dA, dB);
+        require(isDimsCompatibleWithCovers(dA, dB));
+    }
+
+    void init(const Envelope& envA, const Envelope& envB) override {
+        BasicPredicate::requireCovers(envA, envB);
+
+    }
+
+    bool isDetermined() const override {
+        return intersectsExteriorOf(RelateGeometry::GEOM_A);
+    }
+
+    bool valueIM() override {
+        return intMatrix.isCovers();
+    }
+};
+
+static std::unique_ptr<IMPredicate> covers();
+
+
+/************************************************************************
+* Creates a predicate to determine whether a geometry is covered
+* by another geometry.
+*
+* The coveredBy predicate has the following equivalent definitions:
+*
+* Every point of this geometry is a point of the other geometry.
+* The DE-9IM Intersection Matrix for the two geometries matches
+* at least one of the following patterns:
+*
+*   [T*F**F***]
+*   [*TF**F***]
+*   [**FT*F***]
+*   [**F*TF***]
+*
+* covers(B, A) = true
+* (coveredBy is the converse of covers())
+*
+* If either geometry is empty, the value of this predicate is false.
+*
+* This predicate is similar to within(),
+* but is more inclusive (i.e. returns true for more cases).
+*
+* @return the predicate instance
+*
+* @see #covers()
+*/
+class CoveredByPredicate : public IMPredicate {
+
+    std::string name() const override {
+        return std::string("coveredBy");
+    }
+
+    bool requireCovers(bool isSourceA) override {
+        return isSourceA == RelateGeometry::GEOM_B;
+    }
+
+    bool requireExteriorCheck(bool isSourceA) const override {
+        //-- only need to check B against Exterior of A
+        return isSourceA == RelateGeometry::GEOM_A;
+    }
+
+    void init(int dA, int dB) override {
+        IMPredicate::init(dA, dB);
+        require(isDimsCompatibleWithCovers(dB, dA));
+    }
+
+    void init(const Envelope& envA, const Envelope& envB) override {
+        BasicPredicate::requireCovers(envB, envA);
+    }
+
+    bool isDetermined() const override {
+        return intersectsExteriorOf(RelateGeometry::GEOM_B);
+    }
+
+    bool valueIM() override {
+        return intMatrix.isCoveredBy();
+    }
+
+};
+
+static std::unique_ptr<IMPredicate> coveredBy();
+
+
+/************************************************************************
+* Creates a predicate to determine whether a geometry crosses another geometry.
+*
+* The crosses predicate has the following equivalent definitions:
+*
+* The geometries have some but not all interior points in common.
+* The DE-9IM Intersection Matrix for the two geometries matches
+* one of the following patterns:
+*
+*    [T*T******] (for P/L, P/A, and L/A cases)
+*    [T*****T**] (for L/P, A/P, and A/L cases)
+*    [0********] (for L/L cases)
+*
+*
+* For the A/A and P/P cases this predicate returns false.
+*
+* The SFS defined this predicate only for P/L, P/A, L/L, and L/A cases.
+* To make the relation symmetric
+* JTS extends the definition to apply to L/P, A/P and A/L cases as well.
+*
+* @return the predicate instance
+*/
+
+class CrossesPredicate : public IMPredicate {
+
+    std::string name() const override {
+        return std::string("crosses");
+    }
+
+    void init(int dA, int dB) override {
+        IMPredicate::init(dA, dB);
+        bool isBothPointsOrAreas =
+            (dA == Dimension::P && dB == Dimension::P) ||
+            (dA == Dimension::A && dB == Dimension::A);
+        require(!isBothPointsOrAreas);
+    }
+
+    bool isDetermined() const override {
+        if (dimA == Dimension::L && dimB == Dimension::L) {
+            //-- L/L interaction can only be dim = P
+            if (getDimension(Location::INTERIOR, Location::INTERIOR) > Dimension::P)
+                return true;
+        }
+        else if (dimA < dimB) {
+            if (isIntersects(Location::INTERIOR, Location::INTERIOR) &&
+                isIntersects(Location::INTERIOR, Location::EXTERIOR)) {
+                return true;
+            }
+        }
+        else if (dimA > dimB) {
+            if (isIntersects(Location::INTERIOR, Location::INTERIOR) &&
+                isIntersects(Location::EXTERIOR, Location::INTERIOR)) {
+                return true;
+            }
+        }
+        return false;
+    }
+
+    bool valueIM() override {
+        return intMatrix.isCrosses(dimA, dimB);
+    }
+};
+
+static std::unique_ptr<IMPredicate> crosses();
+
+
+/************************************************************************
+* Creates a predicate to determine whether two geometries are
+* topologically equal.
+*
+* The equals predicate has the following equivalent definitions:
+*
+* The two geometries have at least one point in common,
+* and no point of either geometry lies in the exterior of the other geometry.
+* The DE-9IM Intersection Matrix for the two geometries matches
+* the pattern T*F**FFF*
+*
+* @return the predicate instance
+*/
+class EqualsTopoPredicate : public IMPredicate {
+
+    std::string name() const override {
+        return std::string("equals");
+    }
+
+    void init(int dA, int dB) override {
+        IMPredicate::init(dA, dB);
+        require(dA == dB);
+    }
+
+    void init(const Envelope& envA, const Envelope& envB) override {
+        require(envA.equals(&envB));
+    }
+
+    bool isDetermined() const override {
+        bool isEitherExteriorIntersects =
+            isIntersects(Location::INTERIOR, Location::EXTERIOR) ||
+            isIntersects(Location::BOUNDARY, Location::EXTERIOR) ||
+            isIntersects(Location::EXTERIOR, Location::INTERIOR) ||
+            isIntersects(Location::EXTERIOR, Location::BOUNDARY);
+
+        return isEitherExteriorIntersects;
+    }
+
+    bool valueIM() override {
+        return intMatrix.isEquals(dimA, dimB);
+    }
+
+};
+
+static std::unique_ptr<IMPredicate> equalsTopo();
+
+
+/************************************************************************
+ * Creates a predicate to determine whether a geometry overlaps another geometry.
+ *
+ * The overlaps predicate has the following equivalent definitions:
+ *
+ * The geometries have at least one point each not shared by the other
+ *     (or equivalently neither covers the other),
+ *     they have the same dimension,
+ *     and the intersection of the interiors of the two geometries has
+ *     the same dimension as the geometries themselves.
+ * The DE-9IM Intersection Matrix for the two geometries matches
+ *     [T*T***T**] (for P/P and A/A cases)
+ *     or [1*T***T**] (for L/L cases)
+ *
+ * If the geometries are of different dimension this predicate returns false.
+ * This predicate is symmetric.
+ *
+ * @return the predicate instance
+ */
+class OverlapsPredicate : public IMPredicate {
+
+    std::string name() const override {
+        return std::string("overlaps");
+    }
+
+    void init(int dA, int dB) override {
+        IMPredicate::init(dA, dB);
+        require(dA == dB);
+    }
+
+    bool isDetermined() const override {
+        if (dimA == Dimension::A || dimA == Dimension::P) {
+            if (isIntersects(Location::INTERIOR, Location::INTERIOR) &&
+                isIntersects(Location::INTERIOR, Location::EXTERIOR) &&
+                isIntersects(Location::EXTERIOR, Location::INTERIOR))
+            return true;
+        }
+        if (dimA == Dimension::L) {
+            if (isDimension(Location::INTERIOR, Location::INTERIOR, Dimension::L) &&
+                isIntersects(Location::INTERIOR, Location::EXTERIOR) &&
+                isIntersects(Location::EXTERIOR, Location::INTERIOR))
+            return true;
+        }
+        return false;
+    }
+
+    bool valueIM() override {
+        return intMatrix.isOverlaps(dimA, dimB);
+    }
+};
+
+static std::unique_ptr<IMPredicate> overlaps();
+
+
+
+
+
+/************************************************************************
+* Creates a predicate to determine whether a geometry touches another geometry.
+*
+* The touches predicate has the following equivalent definitions:
+*
+* The geometries have at least one point in common,
+* but their interiors do not intersect.
+* The DE-9IM Intersection Matrix for the two geometries matches
+* at least one of the following patterns
+*
+*   [FT*******]
+*   [F**T*****]
+*   [F***T****]
+*
+*
+* If both geometries have dimension 0, the predicate returns false,
+* since points have only interiors.
+* This predicate is symmetric.
+*
+* @return the predicate instance
+*/
+class TouchesPredicate : public IMPredicate {
+
+    std::string name() const override {
+        return std::string("touches");
+    }
+
+    void init(int dA, int dB) override {
+        IMPredicate::init(dA, dB);
+        bool isBothPoints = (dA == 0 && dB == 0);
+        require(! isBothPoints);
+    }
+
+    bool isDetermined() const override {
+        bool isInteriorsIntersects = isIntersects(Location::INTERIOR, Location::INTERIOR);
+        return isInteriorsIntersects;
+    }
+
+    bool valueIM() override {
+        return intMatrix.isTouches(dimA, dimB);
+    }
+};
+
+static std::unique_ptr<IMPredicate> touches();
+
+/**
+ * Creates a predicate that matches a DE-9IM matrix pattern.
+ *
+ * @param imPattern the pattern to match
+ * @return a predicate that matches the pattern
+ *
+ * @see IntersectionMatrixPattern
+ */
+static std::unique_ptr<TopologyPredicate> matches(const std::string& imPattern)
+{
+    return std::unique_ptr<TopologyPredicate>(new IMPatternMatcher(imPattern));
+}
+
+
+
+}; // !RelatePredicate
+
+
+} // 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
new file mode 100644
index 000000000..f6f284e40
--- /dev/null
+++ b/include/geos/operation/relateng/RelateSegmentString.h
@@ -0,0 +1,170 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (c) 2024 Martin Davis
+ * Copyright (C) 2024 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#pragma once
+
+#include <geos/noding/BasicSegmentString.h>
+#include <geos/export.h>
+
+
+#include <string>
+#include <sstream>
+
+
+// Forward declarations
+namespace geos {
+namespace geom {
+    class CoordinateXY;
+    class CoordinateSequence;
+    class Geometry;
+}
+namespace operation {
+namespace relateng {
+    class RelateGeometry;
+    class NodeSection;
+}
+}
+}
+
+
+using geos::noding::BasicSegmentString;
+using geos::geom::Geometry;
+using geos::geom::CoordinateXY;
+using geos::geom::CoordinateSequence;
+
+
+namespace geos {      // geos.
+namespace operation { // geos.operation
+namespace relateng { // geos.operation.relateng
+
+
+
+class GEOS_DLL RelateSegmentString : public BasicSegmentString {
+
+private:
+
+    // Members
+    bool m_isA;
+    int m_dimension;
+    int m_id;
+    int m_ringId;
+    const RelateGeometry* m_inputGeom;
+    const Geometry* m_parentPolygonal = nullptr;
+
+    std::unique_ptr<CoordinateSequence> csStore;
+
+    // Constructor
+    RelateSegmentString(
+        const CoordinateSequence* pts,
+        bool isA,
+        int dimension,
+        int id,
+        int ringId,
+        const Geometry* poly,
+        const RelateGeometry* inputGeom,
+        bool orient)
+        : BasicSegmentString(const_cast<CoordinateSequence*>(pts), nullptr)
+        , m_isA(isA)
+        , m_dimension(dimension)
+        , m_id(id)
+        , 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(
+        const CoordinateSequence* pts,
+        bool isA, int dim, int elementId, int ringId,
+        const Geometry* poly, const RelateGeometry* parent, bool orient);
+
+    /**
+     *
+     * @param ss
+     * @param segIndex
+     * @param pt
+     * @return the previous vertex, or null if none exists
+     */
+    const CoordinateXY* prevVertex(
+        std::size_t segIndex,
+        const CoordinateXY* pt) const;
+
+    /**
+     * @param ss
+     * @param segIndex
+     * @param pt
+     * @return the next vertex, or null if none exists
+     */
+    const CoordinateXY* nextVertex(
+        std::size_t segIndex,
+        const CoordinateXY* pt) const;
+
+    void orientAndRemoveRepeated(bool orientCW);
+
+    void removeRepeated();
+
+
+public:
+
+    static std::unique_ptr<RelateSegmentString>
+    createLine(
+        const CoordinateSequence* pts,
+        bool isA, int elementId,
+        const RelateGeometry* parent, bool orient = false);
+
+    static std::unique_ptr<RelateSegmentString>
+    createRing(
+        const CoordinateSequence* pts,
+        bool isA, int elementId, int ringId,
+        const Geometry* poly, const RelateGeometry* parent, bool orient = true);
+
+    bool isA() const;
+
+    const RelateGeometry* getGeometry() const;
+
+    const Geometry* getPolygonal() const;
+
+    NodeSection* createNodeSection(std::size_t segIndex, const CoordinateXY& intPt) const;
+
+    /**
+     * Tests if a segment intersection point has that segment as its
+     * canonical containing segment.
+     * Segments are half-closed, and contain their start point but not the endpoint,
+     * except for the final segment in a non-closed segment string, which contains
+     * its endpoint as well.
+     * This test ensures that vertices are assigned to a unique segment in a segment string.
+     * In particular, this avoids double-counting intersections which lie exactly
+     * at segment endpoints.
+     *
+     * @param segIndex the segment the point may lie on
+     * @param pt the point
+     * @return true if the segment contains the point
+     */
+    bool isContainingSegment(std::size_t segIndex, const CoordinateXY* pt) const;
+
+};
+
+} // namespace geos.operation.relateng
+} // namespace geos.operation
+} // namespace geos
+
diff --git a/include/geos/operation/relateng/TopologyPredicate.h b/include/geos/operation/relateng/TopologyPredicate.h
new file mode 100644
index 000000000..ce6e27e47
--- /dev/null
+++ b/include/geos/operation/relateng/TopologyPredicate.h
@@ -0,0 +1,199 @@
+/**********************************************************************
+ *
+ * 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/geom/Location.h>
+#include <geos/export.h>
+
+
+// Forward declarations
+namespace geos {
+namespace geom {
+    class Envelope;
+}
+}
+
+
+using geos::geom::Envelope;
+using geos::geom::Location;
+
+
+namespace geos {      // geos.
+namespace operation { // geos.operation
+namespace relateng {  // geos.operation.relateng
+
+
+class GEOS_DLL TopologyPredicate {
+
+public:
+
+    /* Virtual destructor to ensure proper cleanup of derived classes */
+    virtual ~TopologyPredicate() {};
+
+    /**
+     * Gets the name of the predicate.
+     *
+     * @return the predicate name
+     */
+    virtual std::string name() const = 0;
+
+    /**
+     * Indicates that the value of the predicate can be finalized
+     * based on its current state.
+     */
+    virtual void finish() = 0;
+
+    /**
+     * Tests if the predicate value is known.
+     *
+     * @return true if the result is known
+     */
+    virtual bool isKnown() const = 0;
+
+    /**
+     * Gets the current value of the predicate result.
+     * The value is only valid if isKnown() is true.
+     *
+     * @return the predicate result value
+     */
+    virtual bool value() const = 0;
+
+    /**
+     * Reports whether this predicate requires self-noding for
+     * geometries which contain crossing edges
+     * (for example, LineString, or GeometryCollection
+     * containing lines or polygons which may self-intersect).
+     * Self-noding ensures that intersections are computed consistently
+     * in cases which contain self-crossings and mutual crossings.
+     *
+     * Most predicates require this, but it can
+     * be avoided for simple intersection detection
+     * (such as in RelatePredicate#intersects()
+     * and RelatePredicate#disjoint().
+     * Avoiding self-noding improves performance for polygonal inputs.
+     *
+     * @return true if self-noding is required.
+     */
+    virtual bool requireSelfNoding() const {
+        return true;
+    };
+
+    /**
+     * Reports whether this predicate requires interaction between
+     * the input geometries.
+     * This is the case if
+     *
+     * IM[I, I] >= 0 or IM[I, B] >= 0 or IM[B, I] >= 0 or IM[B, B] >= 0
+     *
+     * This allows a fast result if
+     * the envelopes of the geometries are disjoint.
+     *
+     * @return true if the geometries must interact
+     */
+    virtual bool requireInteraction() const {
+        return true;
+    };
+
+    /**
+     * Reports whether this predicate requires that the source
+     * cover the target.
+     * This is the case if
+     *
+     * IM[Ext(Src), Int(Tgt)] = F and IM[Ext(Src), Bdy(Tgt)] = F
+     *
+     * If true, this allows a fast result if
+     * the source envelope does not cover the target envelope.
+     *
+     * @param isSourceA indicates the source input geometry
+     * @return true if the predicate requires checking whether the source covers the target
+     */
+    virtual bool requireCovers(bool isSourceA) {
+        (void)isSourceA;
+        return false;
+    }
+
+    /**
+     * Reports whether this predicate requires checking if the source input intersects
+     * the Exterior of the target input.
+     * This is the case if:
+     *
+     * IM[Int(Src), Ext(Tgt)] >= 0 or IM[Bdy(Src), Ext(Tgt)] >= 0
+     *
+     * If false, this may permit a faster result in some geometric situations.
+     *
+     * @param isSourceA indicates the source input geometry
+     * @return true if the predicate requires checking whether the source intersects the target exterior
+     */
+    virtual bool requireExteriorCheck(bool isSourceA) const {
+        (void)isSourceA;
+        return true;
+    }
+
+    /**
+     * Initializes the predicate for a specific geometric case.
+     * This may allow the predicate result to become known
+     * if it can be inferred from the dimensions.
+     *
+     * @param dimA the dimension of geometry A
+     * @param dimB the dimension of geometry B
+     *
+     * @see Dimension
+     */
+    virtual void init(int dimA, int dimB) {
+        (void)dimA;
+        (void)dimB;
+    };
+
+    /**
+     * Initializes the predicate for a specific geometric case.
+     * This may allow the predicate result to become known
+     * if it can be inferred from the envelopes.
+     *
+     * @param envA the envelope of geometry A
+     * @param envB the envelope of geometry B
+     */
+    virtual void init(const Envelope& envA, const Envelope& envB)
+    {
+        //-- default if envelopes provide no information
+        (void)envA;
+        (void)envB;
+    };
+
+    /**
+     * Updates the entry in the DE-9IM intersection matrix
+     * for given Location in the input geometries.
+     *
+     * If this method is called with a {@link Dimension} value
+     * which is less than the current value for the matrix entry,
+     * the implementing class should avoid changing the entry
+     * if this would cause information loss.
+     *
+     * @param locA the location on the A axis of the matrix
+     * @param locB the location on the B axis of the matrix
+     * @param dimension the dimension value for the entry
+     *
+     * @see Dimension
+     * @see Location
+     */
+    virtual void updateDimension(Location locA, Location locB, int dimension) = 0;
+
+
+};
+
+} // namespace geos.operation.relateng
+} // namespace geos.operation
+} // namespace geos
+
diff --git a/src/geom/GeometryCollection.cpp b/src/geom/GeometryCollection.cpp
index 096ca3526..b1497eb4b 100644
--- a/src/geom/GeometryCollection.cpp
+++ b/src/geom/GeometryCollection.cpp
@@ -456,5 +456,24 @@ 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
 } // namespace geos
diff --git a/src/operation/relateng/BasicPredicate.cpp b/src/operation/relateng/BasicPredicate.cpp
new file mode 100644
index 000000000..c225de6fe
--- /dev/null
+++ b/src/operation/relateng/BasicPredicate.cpp
@@ -0,0 +1,137 @@
+/**********************************************************************
+ *
+ * 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/operation/relateng/BasicPredicate.h>
+#include <geos/geom/Envelope.h>
+#include <geos/geom/Location.h>
+#include <geos/constants.h>
+
+#include <sstream>
+
+
+using geos::geom::Envelope;
+using geos::geom::Location;
+
+
+namespace geos {      // geos
+namespace operation { // geos.operation
+namespace relateng {  // geos.operation.relateng
+
+
+
+/* private static */
+bool
+BasicPredicate::isKnown(int val)
+{
+     return val > UNKNOWN;
+}
+
+/* private static */
+bool
+BasicPredicate::toBoolean(int val)
+{
+    return val == TRUE;
+}
+
+/* private static */
+int
+BasicPredicate::toValue(bool val)
+{
+    return val ? TRUE : FALSE;
+}
+
+
+/* public static */
+bool
+BasicPredicate::isIntersection(Location locA, Location locB)
+{
+    //-- i.e. some location on both geometries intersects
+    return locA != Location::EXTERIOR && locB != Location::EXTERIOR;
+}
+
+
+// /* public */
+// bool isSelfNodingRequired() {
+//     return false;
+//   }
+
+
+/* public override */
+bool
+BasicPredicate::isKnown() const
+{
+    return isKnown(m_value);
+}
+
+/* public override */
+bool
+BasicPredicate::value() const
+{
+    return toBoolean(m_value);
+}
+
+
+/* protected */
+void
+BasicPredicate::setValue(bool val)
+{
+    //-- don't change already-known value
+    if (isKnown())
+        return;
+    m_value = toValue(val);
+}
+
+/* protected */
+void
+BasicPredicate::setValue(int val)
+{
+    //-- don't change already-known value
+    if (isKnown())
+        return;
+    m_value = val;
+}
+
+
+/* protected */
+void
+BasicPredicate::setValueIf(bool val, bool cond)
+{
+    if (cond)
+        setValue(val);
+}
+
+/* protected */
+void
+BasicPredicate::require(bool cond)
+{
+    if (! cond)
+        setValue(false);
+}
+
+/* protected */
+void
+BasicPredicate::requireCovers(const Envelope& a, const Envelope& b)
+{
+    require(a.covers(b));
+}
+
+
+} // namespace geos.operation.overlayng
+} // namespace geos.operation
+} // namespace geos
+
+
+
+
diff --git a/src/operation/relateng/IMPatternMatcher.cpp b/src/operation/relateng/IMPatternMatcher.cpp
new file mode 100644
index 000000000..d70c88044
--- /dev/null
+++ b/src/operation/relateng/IMPatternMatcher.cpp
@@ -0,0 +1,149 @@
+/**********************************************************************
+ *
+ * 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/operation/relateng/IMPatternMatcher.h>
+#include <geos/geom/Envelope.h>
+#include <geos/geom/Location.h>
+
+#include <sstream>
+
+
+using geos::geom::Envelope;
+using geos::geom::Location;
+
+
+namespace geos {      // geos
+namespace operation { // geos.operation
+namespace relateng {  // geos.operation.relateng
+
+
+/* public */
+std::string
+IMPatternMatcher::name() const
+{
+    return "IMPattern";
+}
+
+
+/* public */
+void
+IMPatternMatcher::init(const Envelope& envA, const Envelope& envB)
+{
+    IMPredicate::init(dimA, dimB);
+    //-- if pattern specifies any non-E/non-E interaction, envelopes must not be disjoint
+    bool requiresInteraction = requireInteraction(patternMatrix);
+    bool isDisjoint = envA.disjoint(&envB);
+    setValueIf(false, requiresInteraction && isDisjoint);
+}
+
+
+/* public */
+bool
+IMPatternMatcher::requireInteraction() const
+{
+    return requireInteraction(patternMatrix);
+}
+
+
+/* private static */
+bool
+IMPatternMatcher::requireInteraction(const IntersectionMatrix& im)
+{
+    bool requiresInteraction =
+        isInteraction(im.get(Location::INTERIOR, Location::INTERIOR)) ||
+        isInteraction(im.get(Location::INTERIOR, Location::BOUNDARY)) ||
+        isInteraction(im.get(Location::BOUNDARY, Location::INTERIOR)) ||
+        isInteraction(im.get(Location::BOUNDARY, Location::BOUNDARY));
+    return requiresInteraction;
+}
+
+
+/* private static */
+bool
+IMPatternMatcher::isInteraction(int imDim)
+{
+    return imDim == Dimension::True || imDim >= Dimension::P;
+}
+
+
+/* public */
+bool
+IMPatternMatcher::isDetermined() const
+{
+    /**
+     * Matrix entries only increase in dimension as topology is computed.
+     * The predicate can be short-circuited (as false) if
+     * any computed entry is greater than the mask value.
+     */
+    std::array<Location,3> locs = {
+        Location::INTERIOR, Location::BOUNDARY, Location::EXTERIOR};
+
+    for (Location i : locs) {
+        for (Location j : locs) {
+            int patternEntry = patternMatrix.get(i, j);
+
+            if (patternEntry == Dimension::DONTCARE)
+                continue;
+
+            int matrixVal = getDimension(i, j);
+
+            //-- mask entry TRUE requires a known matrix entry
+            if (patternEntry == Dimension::True) {
+                if (matrixVal < 0)
+                    return false;
+            }
+            //-- result is known (false) if matrix entry has exceeded mask
+            else if (matrixVal > patternEntry)
+                return true;
+        }
+    }
+    return false;
+}
+
+
+/* public */
+bool
+IMPatternMatcher::valueIM()
+{
+    bool val = intMatrix.matches(imPattern);
+    return val;
+}
+
+
+/* public */
+std::string
+IMPatternMatcher::toString() const
+{
+    return name() + "(" + imPattern + ")";
+}
+
+
+/* public friend */
+std::ostream&
+operator<<(std::ostream& os, const IMPatternMatcher& imp)
+{
+    os << imp.toString();
+    return os;
+}
+
+
+
+} // namespace geos.operation.overlayng
+} // namespace geos.operation
+} // namespace geos
+
+
+
+
diff --git a/src/operation/relateng/IMPredicate.cpp b/src/operation/relateng/IMPredicate.cpp
new file mode 100644
index 000000000..c41be5e8c
--- /dev/null
+++ b/src/operation/relateng/IMPredicate.cpp
@@ -0,0 +1,155 @@
+/**********************************************************************
+ *
+ * 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/operation/relateng/IMPredicate.h>
+#include <geos/geom/Envelope.h>
+#include <geos/geom/Location.h>
+#include <geos/constants.h>
+
+#include <sstream>
+
+
+using geos::geom::Envelope;
+using geos::geom::Location;
+
+
+namespace geos {      // geos
+namespace operation { // geos.operation
+namespace relateng {  // geos.operation.relateng
+
+
+/* public static */
+bool
+IMPredicate::isDimsCompatibleWithCovers(int dim0, int dim1)
+{
+    //- allow Points coveredBy zero-length Lines
+    if (dim0 == Dimension::P && dim1 == Dimension::L)
+        return true;
+    return dim0 >= dim1;
+}
+
+
+/* public */
+void
+IMPredicate::init(int dA, int dB)
+{
+    dimA = dA;
+    dimB = dB;
+}
+
+
+/* public */
+void
+IMPredicate::updateDimension(Location locA, Location locB, int dimension)
+{
+    //-- only record an increased dimension value
+    if (isDimChanged(locA, locB, dimension)) {
+        intMatrix.set(locA, locB, dimension);
+        //-- set value if predicate value can be known
+        if (isDetermined()) {
+            setValue(valueIM());
+        }
+    }
+}
+
+
+/* public */
+bool
+IMPredicate::isDimChanged(Location locA, Location locB, int dimension) const
+{
+    return dimension > intMatrix.get(locA, locB);
+}
+
+
+/* protected */
+bool
+IMPredicate::intersectsExteriorOf(bool isA) const
+{
+    if (isA) {
+        return isIntersects(Location::EXTERIOR, Location::INTERIOR)
+            || isIntersects(Location::EXTERIOR, Location::BOUNDARY);
+    }
+    else {
+        return isIntersects(Location::INTERIOR, Location::EXTERIOR)
+            || isIntersects(Location::BOUNDARY, Location::EXTERIOR);
+    }
+}
+
+
+/* protected */
+bool
+IMPredicate::isIntersects(Location locA, Location locB) const
+{
+    return intMatrix.get(locA, locB) >= Dimension::P;
+}
+
+
+/* public */
+bool
+IMPredicate::isKnown(Location locA, Location locB) const
+{
+    return intMatrix.get(locA, locB) != DIM_UNKNOWN;
+}
+
+
+/* public */
+bool
+IMPredicate::isDimension(Location locA, Location locB, int dimension) const
+{
+    return intMatrix.get(locA, locB) == dimension;
+}
+
+
+/* public */
+int
+IMPredicate::getDimension(Location locA, Location locB) const
+{
+    return intMatrix.get(locA, locB);
+}
+
+
+/* public */
+void
+IMPredicate::finish()
+{
+    setValue(valueIM());
+}
+
+
+/* public */
+std::string
+IMPredicate::toString() const
+{
+    return name() + ": " + intMatrix.toString();
+}
+
+
+/* public friend */
+std::ostream&
+operator<<(std::ostream& os, const IMPredicate& imp)
+{
+    os << imp.toString();
+    return os;
+}
+
+
+
+} // namespace geos.operation.overlayng
+} // namespace geos.operation
+} // namespace geos
+
+
+
+
diff --git a/src/operation/relateng/LinearBoundary.cpp b/src/operation/relateng/LinearBoundary.cpp
index 5b4dfbae7..c97d6b84b 100644
--- a/src/operation/relateng/LinearBoundary.cpp
+++ b/src/operation/relateng/LinearBoundary.cpp
@@ -100,9 +100,6 @@ LinearBoundary::addEndpoint(const CoordinateXY *p, Coordinate::ConstIntMap& vert
     if (it != vertexDegree.end()) {
         dim = it->second;
     }
-    // dim++;
-    // std::pair<const CoordinateXY*, int> entry(p, dim);
-    // vertexDegree.insert(entry);
     vertexDegree[p] = dim + 1;
 }
 
diff --git a/src/operation/relateng/NodeSection.cpp b/src/operation/relateng/NodeSection.cpp
index 3fdfac851..81d7d304e 100644
--- a/src/operation/relateng/NodeSection.cpp
+++ b/src/operation/relateng/NodeSection.cpp
@@ -20,9 +20,6 @@
 #include <geos/operation/relateng/NodeSection.h>
 #include <sstream>
 
-namespace geos {      // geos
-namespace operation { // geos.operation
-namespace relateng {  // geos.operation.relateng
 
 using geos::geom::Coordinate;
 using geos::geom::CoordinateXY;
@@ -31,6 +28,11 @@ using geos::geom::Dimension;
 using geos::io::WKTWriter;
 
 
+namespace geos {      // geos
+namespace operation { // geos.operation
+namespace relateng {  // geos.operation.relateng
+
+
 /* public */
 const CoordinateXY *
 NodeSection::getVertex(int i) const
diff --git a/src/operation/relateng/PolygonNodeConverter.cpp b/src/operation/relateng/PolygonNodeConverter.cpp
index 54c080956..e9b1e9c41 100644
--- a/src/operation/relateng/PolygonNodeConverter.cpp
+++ b/src/operation/relateng/PolygonNodeConverter.cpp
@@ -43,10 +43,11 @@ PolygonNodeConverter::convert(std::vector<const NodeSection*>& polySections)
         const NodeSection* ns1,
         const NodeSection* ns2)
     {
-        return PolygonNodeTopology::compareAngle(
+        int comp = PolygonNodeTopology::compareAngle(
             ns1->nodePt(),
             ns1->getVertex(0),
             ns2->getVertex(0));
+        return comp < 0;
     };
 
     std::sort(polySections.begin(), polySections.end(), comparator);
diff --git a/src/operation/relateng/RelateEdge.cpp b/src/operation/relateng/RelateEdge.cpp
index b92cd93c3..4670b058b 100644
--- a/src/operation/relateng/RelateEdge.cpp
+++ b/src/operation/relateng/RelateEdge.cpp
@@ -27,7 +27,6 @@
 using geos::algorithm::PolygonNodeTopology;
 using geos::geom::CoordinateXY;
 using geos::geom::Dimension;
-using geos::geom::Location;
 using geos::geom::Position;
 using geos::io::WKTWriter;
 
diff --git a/src/operation/relateng/RelateGeometry.cpp b/src/operation/relateng/RelateGeometry.cpp
new file mode 100644
index 000000000..d31acf30e
--- /dev/null
+++ b/src/operation/relateng/RelateGeometry.cpp
@@ -0,0 +1,464 @@
+/**********************************************************************
+ *
+ * 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/GeometryCollection.h>
+#include <geos/geom/LineString.h>
+#include <geos/geom/LinearRing.h>
+#include <geos/geom/Location.h>
+#include <geos/geom/MultiPolygon.h>
+#include <geos/geom/Point.h>
+#include <geos/geom/Polygon.h>
+#include <geos/geom/util/ComponentCoordinateExtracter.h>
+#include <geos/geom/util/PointExtracter.h>
+#include <geos/operation/relateng/RelateGeometry.h>
+#include <geos/operation/relateng/RelateSegmentString.h>
+#include <geos/operation/relateng/DimensionLocation.h>
+
+
+#include <sstream>
+
+
+using geos::algorithm::BoundaryNodeRule;
+using namespace geos::geom;
+
+
+namespace geos {      // geos
+namespace operation { // geos.operation
+namespace relateng {  // geos.operation.relateng
+
+
+RelateGeometry::RelateGeometry(const Geometry* input, bool isPrepared, const BoundaryNodeRule& bnRule)
+    : geom(input)
+    , m_isPrepared(isPrepared)
+    , geomEnv(input->getEnvelopeInternal())
+    , boundaryNodeRule(bnRule)
+    , geomDim(input->getDimension())
+    , isLineZeroLen(isZeroLength(input))
+    , isGeomEmpty(input->isEmpty())
+{
+    analyzeDimensions();
+}
+
+
+/* public static */
+std::string
+RelateGeometry::name(bool isA)
+{
+    return isA ? "A" : "B";
+}
+
+
+/* private */
+void
+RelateGeometry::analyzeDimensions()
+{
+    if (isGeomEmpty) {
+        return;
+    }
+    GeometryTypeId typeId = geom->getGeometryTypeId();
+    if (typeId == GEOS_POINT || typeId == GEOS_MULTIPOINT) {
+        hasPoints = true;
+        geomDim = Dimension::P;
+        return;
+    }
+    if (typeId == GEOS_LINESTRING || typeId == GEOS_MULTILINESTRING) {
+        hasLines = true;
+        geomDim = Dimension::L;
+        return;
+    }
+    if (typeId == GEOS_POLYGON || typeId == GEOS_MULTIPOLYGON) {
+        hasAreas = true;
+        geomDim = Dimension::A;
+        return;
+    }
+    //-- analyze a (possibly mixed type) collection
+    std::vector<const Geometry*> elems;
+    const GeometryCollection* col = static_cast<const GeometryCollection*>(geom);
+    col->getAllGeometries(elems);
+    for (const Geometry* elem : elems)
+    {
+        if (elem->isEmpty())
+            continue;
+        if (elem->getGeometryTypeId() == GEOS_POINT) {
+            hasPoints = true;
+            if (geomDim < Dimension::P) geomDim = Dimension::P;
+        }
+        if (elem->getGeometryTypeId() == GEOS_LINESTRING) {
+            hasLines = true;
+            if (geomDim < Dimension::L) geomDim = Dimension::L;
+        }
+        if (elem->getGeometryTypeId() == GEOS_POLYGON) {
+            hasAreas = true;
+            if (geomDim < Dimension::A) geomDim = Dimension::A;
+        }
+    }
+}
+
+
+/* private static */
+bool
+RelateGeometry::isZeroLength(const Geometry* geom)
+{
+    std::vector<const Geometry*> elems;
+    const GeometryCollection* col = static_cast<const GeometryCollection*>(geom);
+    col->getAllGeometries(elems);
+    for (const Geometry* elem : elems) {
+        if (elem->getGeometryTypeId() == GEOS_LINESTRING) {
+            if (! isZeroLength(static_cast<const LineString*>(elem))) {
+                return false;
+            }
+        }
+    }
+    return true;
+}
+
+/* private static */
+bool
+RelateGeometry::isZeroLength(const LineString* line) {
+    if (line->getNumPoints() >= 2) {
+        const CoordinateXY& p0 = line->getCoordinateN(0);
+        for (std::size_t i = 1; i < line->getNumPoints(); i++) {
+            // NOTE !!! CHANGE FROM JTS, original below
+            // const CoordinateXY& pi = line.getCoordinateN(1);
+            const CoordinateXY& pi = line->getCoordinateN(i);
+            //-- most non-zero-len lines will trigger this right away
+            if (! p0.equals2D(pi)) {
+                return false;
+            }
+        }
+    }
+    return true;
+}
+
+
+/* public */
+const Geometry*
+RelateGeometry::getGeometry() const
+{
+    return geom;
+}
+
+/* public */
+bool
+RelateGeometry::isPrepared() const
+{
+    return m_isPrepared;
+}
+
+/* public */
+const Envelope*
+RelateGeometry::getEnvelope() const
+{
+    return geomEnv;
+}
+
+/* public */
+int
+RelateGeometry::getDimension() const
+{
+    return geomDim;
+}
+
+/* public */
+bool
+RelateGeometry::hasDimension(int dim) const
+{
+    switch (dim) {
+        case Dimension::P: return hasPoints;
+        case Dimension::L: return hasLines;
+        case Dimension::A: return hasAreas;
+    }
+    return false;
+}
+
+
+/* public */
+int
+RelateGeometry::getDimensionReal() const
+{
+    if (isGeomEmpty)
+        return Dimension::False;
+    if (getDimension() == Dimension::L && isLineZeroLen)
+        return Dimension::P;
+    if (hasAreas)
+        return Dimension::A;
+    if (hasLines)
+        return Dimension::L;
+    return Dimension::P;
+}
+
+/* public */
+bool
+RelateGeometry::hasEdges() const
+{
+    return hasLines || hasAreas;
+}
+
+/* private */
+RelatePointLocator*
+RelateGeometry::getLocator()
+{
+    if (locator == nullptr)
+        locator.reset(new RelatePointLocator(geom, m_isPrepared, boundaryNodeRule));
+    return locator.get();
+}
+
+
+/* public */
+bool
+RelateGeometry::isNodeInArea(const CoordinateXY* nodePt, const Geometry* parentPolygonal)
+{
+    int dimLoc = getLocator()->locateNodeWithDim(nodePt, parentPolygonal);
+    return dimLoc == DimensionLocation::AREA_INTERIOR;
+}
+
+
+/* public */
+Location
+RelateGeometry::locateLineEnd(const CoordinateXY* p)
+{
+    return getLocator()->locateLineEnd(p);
+}
+
+
+/* public */
+Location
+RelateGeometry::locateAreaVertex(const CoordinateXY* pt)
+{
+    /**
+     * Can pass a null polygon, because the point is an exact vertex,
+     * which will be detected as being on the boundary of its polygon
+     */
+    return locateNode(pt, nullptr);
+}
+
+
+/* public */
+Location
+RelateGeometry::locateNode(const CoordinateXY* pt, const Geometry* parentPolygonal)
+{
+    return getLocator()->locateNode(pt, parentPolygonal);
+}
+
+
+/* public */
+int
+RelateGeometry::locateWithDim(const CoordinateXY* pt)
+{
+    int loc = getLocator()->locateWithDim(pt);
+    return loc;
+}
+
+
+/* public */
+bool
+RelateGeometry::isPointsOrPolygons() const
+{
+    GeometryTypeId typeId = geom->getGeometryTypeId();
+    return typeId == GEOS_POINT
+        || typeId == GEOS_MULTIPOINT
+        || typeId == GEOS_POLYGON
+        || typeId == GEOS_MULTIPOLYGON;
+}
+
+
+/* public */
+bool
+RelateGeometry::isPolygonal() const
+{
+    //TODO: also true for a GC containing one polygonal element (and possibly some lower-dimension elements)
+    GeometryTypeId typeId = geom->getGeometryTypeId();
+    return typeId == GEOS_POLYGON
+        || typeId == GEOS_MULTIPOLYGON;
+}
+
+
+/* public */
+bool
+RelateGeometry::isEmpty() const
+{
+    return isGeomEmpty;
+}
+
+
+/* public */
+bool
+RelateGeometry::hasBoundary()
+{
+    return getLocator()->hasBoundary();
+}
+
+
+/* public */
+Coordinate::ConstXYSet&
+RelateGeometry::getUniquePoints()
+{
+    if (uniquePoints.empty()) {
+        uniquePoints = createUniquePoints();
+    }
+    return uniquePoints;
+}
+
+
+/* private */
+Coordinate::ConstXYSet
+RelateGeometry::createUniquePoints()
+{
+    //-- only called on P geometries
+    std::vector<const CoordinateXY*> pts;
+    geom::util::ComponentCoordinateExtracter::getCoordinates(*geom, pts);
+    Coordinate::ConstXYSet set(pts.begin(), pts.end());
+    return set;
+}
+
+
+/* public */
+std::vector<const Point*>
+RelateGeometry::getEffectivePoints()
+{
+    std::vector<const Point*> ptListAll;
+    geom::util::PointExtracter::getPoints(*geom, ptListAll);
+
+    if (getDimensionReal() <= Dimension::P)
+        return ptListAll;
+
+    //-- only return Points not covered by another element
+    std::vector<const Point*> ptList;
+    for (const Point* p : ptListAll) {
+        int locDim = locateWithDim(p->getCoordinate());
+        if (DimensionLocation::dimension(locDim) == Dimension::P) {
+            ptList.push_back(p);
+        }
+    }
+    return ptList;
+}
+
+
+/* public */
+std::vector<std::unique_ptr<RelateSegmentString>>
+RelateGeometry::extractSegmentStrings(bool isA, const Envelope* env)
+{
+    std::vector<std::unique_ptr<RelateSegmentString>> segStrings;
+    extractSegmentStrings(isA, env, geom, segStrings);
+    return segStrings;
+}
+
+
+/* private */
+void
+RelateGeometry::extractSegmentStrings(bool isA,
+    const Envelope* env, const Geometry* p_geom,
+    std::vector<std::unique_ptr<RelateSegmentString>>& segStrings)
+{
+    //-- record if parent is MultiPolygon
+    const MultiPolygon* parentPolygonal = nullptr;
+    if (p_geom->getGeometryTypeId() == GEOS_MULTIPOLYGON) {
+        parentPolygonal = static_cast<const MultiPolygon*>(p_geom);
+    }
+
+    for (std::size_t i = 0; i < p_geom->getNumGeometries(); i++) {
+        const Geometry* g = p_geom->getGeometryN(i);
+        if (g->getGeometryTypeId() == GEOS_GEOMETRYCOLLECTION) {
+            extractSegmentStrings(isA, env, g, segStrings);
+        }
+        else {
+            extractSegmentStringsFromAtomic(isA, g, parentPolygonal, env, segStrings);
+        }
+    }
+}
+
+
+/* private */
+void
+RelateGeometry::extractSegmentStringsFromAtomic(bool isA,
+    const Geometry* p_geom, const MultiPolygon* parentPolygonal,
+    const Envelope* env,
+    std::vector<std::unique_ptr<RelateSegmentString>>& segStrings)
+{
+    if (p_geom->isEmpty())
+        return;
+
+    bool doExtract = (env == nullptr) || env->intersects(p_geom->getEnvelopeInternal());
+    if (! doExtract)
+        return;
+
+    elementId++;
+    if (p_geom->getGeometryTypeId() == GEOS_LINESTRING) {
+        const LineString* line = static_cast<const LineString*>(p_geom);
+        auto cs = line->getCoordinatesRO();
+        auto ss = RelateSegmentString::createLine(cs, isA, elementId, this);
+        segStrings.emplace_back(ss.release());
+    }
+    else if (p_geom->getGeometryTypeId() == GEOS_POLYGON) {
+        const Polygon* poly = static_cast<const Polygon*>(p_geom);
+        const Geometry* parentPoly;
+        if (parentPolygonal != nullptr)
+            parentPoly = static_cast<const Geometry*>(parentPolygonal);
+        else
+            parentPoly = static_cast<const Geometry*>(poly);
+        extractRingToSegmentString(isA, poly->getExteriorRing(), 0, env, parentPoly, segStrings);
+        for (uint32_t i = 0; i < poly->getNumInteriorRing(); i++) {
+            extractRingToSegmentString(isA, poly->getInteriorRingN(i), static_cast<int>(i+1), env, parentPoly, segStrings);
+        }
+    }
+}
+
+
+/* private */
+void
+RelateGeometry::extractRingToSegmentString(bool isA,
+    const LinearRing* ring, int ringId, const Envelope* env,
+    const Geometry* parentPoly,
+    std::vector<std::unique_ptr<RelateSegmentString>>& 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());
+}
+
+
+/* public */
+std::string
+RelateGeometry::toString() const
+{
+    return geom->toString();
+}
+
+
+/* public friend */
+std::ostream&
+operator<<(std::ostream& os, const RelateGeometry& rg)
+{
+    os << rg.toString();
+    return os;
+}
+
+
+
+} // namespace geos.operation.overlayng
+} // namespace geos.operation
+} // namespace geos
+
+
diff --git a/src/operation/relateng/RelatePredicate.cpp b/src/operation/relateng/RelatePredicate.cpp
new file mode 100644
index 000000000..3efe5537f
--- /dev/null
+++ b/src/operation/relateng/RelatePredicate.cpp
@@ -0,0 +1,109 @@
+/**********************************************************************
+ *
+ * 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/operation/relateng/RelatePredicate.h>
+// #include <geos/constants.h>
+
+#include <memory>
+
+
+// using geos::geom::Envelope;
+// using geos::geom::Location;
+
+
+namespace geos {      // geos
+namespace operation { // geos.operation
+namespace relateng {  // geos.operation.relateng
+
+
+/* public static */
+std::unique_ptr<BasicPredicate>
+RelatePredicate::intersects()
+{
+    return std::unique_ptr<BasicPredicate>(new IntersectsPredicate());
+}
+
+/* public static */
+std::unique_ptr<BasicPredicate>
+RelatePredicate::disjoint()
+{
+    return std::unique_ptr<BasicPredicate>(new DisjointPredicate());
+}
+
+/* public static */
+std::unique_ptr<IMPredicate>
+RelatePredicate::contains()
+{
+    return std::unique_ptr<IMPredicate>(new ContainsPredicate());
+}
+
+/* public static */
+std::unique_ptr<IMPredicate>
+RelatePredicate::within()
+{
+    return std::unique_ptr<IMPredicate>(new WithinPredicate());
+}
+
+/* public static */
+std::unique_ptr<IMPredicate>
+RelatePredicate::covers()
+{
+    return std::unique_ptr<IMPredicate>(new CoversPredicate());
+}
+
+/* public static */
+std::unique_ptr<IMPredicate>
+RelatePredicate::coveredBy()
+{
+    return std::unique_ptr<IMPredicate>(new CoveredByPredicate());
+}
+
+/* public static */
+std::unique_ptr<IMPredicate>
+RelatePredicate::crosses()
+{
+    return std::unique_ptr<IMPredicate>(new CrossesPredicate());
+}
+
+
+/* public static */
+std::unique_ptr<IMPredicate>
+RelatePredicate::equalsTopo()
+{
+    return std::unique_ptr<IMPredicate>(new EqualsTopoPredicate());
+}
+
+
+/* public static */
+std::unique_ptr<IMPredicate>
+RelatePredicate::overlaps()
+{
+    return std::unique_ptr<IMPredicate>(new OverlapsPredicate());
+}
+
+/* public static */
+std::unique_ptr<IMPredicate>
+RelatePredicate::touches()
+{
+    return std::unique_ptr<IMPredicate>(new TouchesPredicate());
+}
+
+} // namespace geos.operation.overlayng
+} // namespace geos.operation
+} // namespace geos
+
+
+
+
diff --git a/src/operation/relateng/RelateSegmentString.cpp b/src/operation/relateng/RelateSegmentString.cpp
new file mode 100644
index 000000000..ccdf46d02
--- /dev/null
+++ b/src/operation/relateng/RelateSegmentString.cpp
@@ -0,0 +1,224 @@
+/**********************************************************************
+ *
+ * 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/Orientation.h>
+#include <geos/geom/Dimension.h>
+#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>
+
+
+using geos::geom::Coordinate;
+using geos::geom::CoordinateXY;
+using geos::geom::Dimension;
+using geos::geom::Geometry;
+using geos::algorithm::Orientation;
+using geos::operation::valid::RepeatedPointRemover;
+
+
+namespace geos {      // geos
+namespace operation { // geos.operation
+namespace relateng {  // geos.operation.relateng
+
+
+/* public static */
+std::unique_ptr<RelateSegmentString>
+RelateSegmentString::createLine(
+    const CoordinateSequence* pts,
+    bool isA, int elementId,
+    const RelateGeometry* parent, bool orient)
+{
+    return createSegmentString(pts, isA, Dimension::L, elementId, -1, nullptr, parent, orient);
+}
+
+
+/* public static */
+std::unique_ptr<RelateSegmentString>
+RelateSegmentString::createRing(
+    const CoordinateSequence* pts,
+    bool isA, int elementId, int ringId,
+    const Geometry* poly, const RelateGeometry* parent, bool orient)
+{
+    return createSegmentString(pts, isA, Dimension::A, elementId, ringId, poly, parent, orient);
+}
+
+
+/* private static */
+std::unique_ptr<RelateSegmentString>
+RelateSegmentString::createSegmentString(
+    const CoordinateSequence* pts,
+    bool isA, int dim, int elementId, int ringId,
+    const Geometry* poly, const RelateGeometry* parent, bool orient)
+{
+    return std::unique_ptr<RelateSegmentString>(new RelateSegmentString(pts, isA, dim, elementId, ringId, poly, parent, orient));
+}
+
+
+/* public */
+bool
+RelateSegmentString::isA() const
+{
+    return m_isA;
+}
+
+
+/* public */
+const RelateGeometry*
+RelateSegmentString::getGeometry() const
+{
+    return m_inputGeom;
+}
+
+
+/* public */
+const Geometry*
+RelateSegmentString::getPolygonal() const
+{
+    return m_parentPolygonal;
+}
+
+
+/* public */
+NodeSection*
+RelateSegmentString::createNodeSection(std::size_t segIndex, const CoordinateXY& intPt) const
+{
+    const CoordinateXY& c0 = getCoordinate(segIndex);
+    const CoordinateXY& c1 = getCoordinate(segIndex + 1);
+    bool isNodeAtVertex = intPt.equals2D(c0) || intPt.equals2D(c1);
+    const CoordinateXY* prev = prevVertex(segIndex, &intPt);
+    const CoordinateXY* next = nextVertex(segIndex, &intPt);
+    NodeSection* a = new NodeSection(m_isA, m_dimension, m_id, m_ringId, m_parentPolygonal, isNodeAtVertex, prev, &intPt, next);
+    return a;
+}
+
+
+/* private */
+const CoordinateXY*
+RelateSegmentString::prevVertex(std::size_t segIndex, const CoordinateXY* pt) const
+{
+    const CoordinateXY& segStart = getCoordinate(segIndex);
+    if (! segStart.equals2D(*pt))
+        return &segStart;
+
+    //-- pt is at segment start, so get previous vertex
+    if (segIndex > 0) {
+        const CoordinateXY& seg = getCoordinate(segIndex - 1);
+        return &seg;
+    }
+
+    if (isClosed())
+        return &(prevInRing(segIndex));
+
+    return nullptr;
+}
+
+
+/* private */
+const CoordinateXY*
+RelateSegmentString::nextVertex(std::size_t segIndex, const CoordinateXY* pt) const
+{
+    const CoordinateXY& segEnd = getCoordinate(segIndex + 1);
+    if (! segEnd.equals2D(*pt))
+        return &segEnd;
+
+    //-- pt is at seg end, so get next vertex
+    if (size() == 2 && segIndex == INDEX_UNKNOWN) {
+        const CoordinateXY& seg = getCoordinate(0);
+        return &seg;
+    }
+
+    if (segIndex < size() - 2) {
+        const CoordinateXY& seg = getCoordinate(segIndex + 2);
+        return &seg;
+    }
+
+    if (isClosed())
+        return &(SegmentString::nextInRing(segIndex + 1));
+
+    //-- segstring is not closed, so there is no next segment
+    return nullptr;
+}
+
+
+/* public */
+bool
+RelateSegmentString::isContainingSegment(std::size_t segIndex, const CoordinateXY* pt) const
+{
+    //-- intersection is at segment start vertex - process it
+    const CoordinateXY& c0 = getCoordinate(segIndex);
+    if (pt->equals2D(c0))
+        return true;
+    const CoordinateXY& c1 = getCoordinate(segIndex+1);
+    if (pt->equals2D(c1)) {
+        bool isFinalSegment = segIndex == size() - 2;
+        if (isClosed() || ! isFinalSegment)
+            return false;
+        //-- for final segment, process intersections with final endpoint
+        return true;
+    }
+    //-- intersection is interior - process it
+    return true;
+}
+
+
+/* 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
+} // namespace geos.operation
+} // namespace geos
+
+
diff --git a/tests/unit/operation/relateng/PolygonNodeConverterTest.cpp b/tests/unit/operation/relateng/PolygonNodeConverterTest.cpp
index 4a4899f70..6037df049 100644
--- a/tests/unit/operation/relateng/PolygonNodeConverterTest.cpp
+++ b/tests/unit/operation/relateng/PolygonNodeConverterTest.cpp
@@ -9,6 +9,7 @@
 #include <geos/io/WKTWriter.h>
 #include <geos/geom/Geometry.h>
 #include <geos/geom/Location.h>
+#include <geos/geom/Dimension.h>
 #include <geos/operation/relateng/PolygonNodeConverter.h>
 #include <geos/algorithm/PolygonNodeTopology.h>
 
@@ -30,67 +31,81 @@ namespace tut {
 // Common data used by all tests
 struct test_polygonnodeconverter_data {
 
-    WKTReader r;
-    // WKTWriter w;
+    NodeSection* sectionShell(double v0x, double v0y, double nx, double ny, double v1x, double v1y) {
+        return section(0, v0x, v0y, nx, ny, v1x, v1y);
+    }
 
-    std::unique_ptr<CoordinateSequence>
-    readPts(const std::string& wkt)
+    NodeSection* sectionHole(double v0x, double v0y, double nx, double ny, double v1x, double v1y) {
+        return section(1, v0x, v0y, nx, ny, v1x, v1y);
+    }
+
+    NodeSection* section(int ringId, double v0x, double v0y, double nx, double ny, double v1x, double v1y) {
+        return new NodeSection(true, Dimension::A, 1, ringId, nullptr, false, 
+            new Coordinate(v0x, v0y), new Coordinate(nx, ny), new Coordinate(v1x, v1y)); 
+    }
+
+    std::vector<const NodeSection*> 
+    toPtrVector(const std::vector<std::unique_ptr<NodeSection>>& input)
     {
-        std::unique_ptr<Geometry> line = r.read(wkt);
-        return line->getCoordinates();
+        std::vector<const NodeSection*> vec;
+        for(std::size_t i = 0, n = input.size(); i < n; ++i) {
+            vec.emplace_back(input[i].get());
+        }
+        return vec;
+    }
+
+    bool checkSectionsEqual(std::vector<const NodeSection*>& ns1, 
+                            std::vector<const NodeSection*>& ns2) {
+        if (ns1.size() != ns2.size())
+            return false;
+        sort(ns1);
+        sort(ns2);
+        for (std::size_t i = 0; i < ns1.size(); i++) {
+            int comp = ns1[i]->compareTo(ns2[i]);
+            if (comp != 0)
+                return false;
+        }
+        return true;
+    }
+    void sort(std::vector<const NodeSection*>& ns) {
+
+        // Comparator lambda for sort support
+        auto comparator = [](
+            const NodeSection* a,
+            const NodeSection* b)
+        {
+            return a->compareTo(b) < 0;
+        };
+            
+        std::sort(ns.begin(), ns.end(), comparator);
+    }
+
+    void 
+    freeNodeSections(std::vector<const NodeSection*>& ns) {
+        for (std::size_t i = 0; i < ns.size(); i++) {
+            delete ns[i]->getVertex(0);
+            delete ns[i]->getVertex(1);
+            delete ns[i]->nodePt();
+            delete ns[i];
+        }
     }
 
     void
-    checkCrossing(const std::string& wktA, const std::string& wktB)
+    checkConversion(std::vector<const NodeSection *>& input, 
+                    std::vector<const NodeSection *>& expected)
     {
-        checkCrossing(wktA, wktB, true);
-    }
-
-    void
-    checkNonCrossing(const std::string& wktA, const std::string& wktB)
-    {
-        checkCrossing(wktA, wktB, false);
-    }
-
-    void
-    checkCrossing(const std::string& wktA, const std::string& wktB, bool isExpected)
-    {
-        auto a = readPts(wktA);
-        auto b = readPts(wktB);
-        // assert: a[1] = b[1]
-        bool isCrossing = PolygonNodeTopology::isCrossing(
-            &(a->getAt<CoordinateXY>(1)),
-            &(a->getAt<CoordinateXY>(0)),
-            &(a->getAt<CoordinateXY>(2)),
-            &(b->getAt<CoordinateXY>(0)),
-            &(b->getAt<CoordinateXY>(2)));
-        ensure("checkCrossing", isCrossing == isExpected);
-    }
-
-    void
-    checkInterior(const std::string& wktA, const std::string& wktB)
-    {
-        checkInteriorSegment(wktA, wktB, true);
-    }
-
-    void
-    checkExterior(const std::string& wktA, const std::string& wktB)
-    {
-        checkInteriorSegment(wktA, wktB, false);
-    }
-
-    void
-    checkInteriorSegment(const std::string& wktA, const std::string& wktB, bool isExpected)
-    {
-        auto a = readPts(wktA);
-        auto b = readPts(wktB);
-        // assert: a[1] = b[1]
-        bool isInterior = PolygonNodeTopology::isInteriorSegment(
-            &(a->getAt<CoordinateXY>(1)),
-            &(a->getAt<CoordinateXY>(0)),
-            &(a->getAt<CoordinateXY>(2)),
-            &(b->getAt<CoordinateXY>(1)));
-        ensure("checkInteriorSegment", isInterior == isExpected);
+        auto actual = PolygonNodeConverter::convert( input );
+        auto actualPtr = toPtrVector(actual);
+        bool isEqual = checkSectionsEqual(actualPtr, expected);
+        freeNodeSections(input);
+        freeNodeSections(expected);
+        /*
+        if (! isEqual) {
+            System.out.println("Expected:" + formatSections(expected));
+            System.out.println("Actual:" + formatSections(actual));      
+        }
+        */
+        ensure("checkConversion", isEqual);
     }
 
 };
@@ -105,68 +120,94 @@ group test_polygonnodeconverter_group("geos::operation::relateng::PolygonNodeCon
 //
 
 
-// testNonCrossing
+// testShells
 template<>
 template<>
 void object::test<1> ()
 {
-    checkCrossing("LINESTRING (500 1000, 1000 1000, 1000 1500)",
-        "LINESTRING (1000 500, 1000 1000, 500 1500)");
+    std::vector<const NodeSection *> input{
+        sectionShell( 1,1, 5,5, 9,9 ),
+        sectionShell( 8,9, 5,5, 6,9 ),
+        sectionShell( 4,9, 5,5, 2,9 )
+    };
+    std::vector<const NodeSection *> expected{
+        sectionShell( 1,1, 5,5, 9,9 ),
+        sectionShell( 8,9, 5,5, 6,9 ),
+        sectionShell( 4,9, 5,5, 2,9 )
+    };
+    checkConversion(input, expected);
 }
 
-// testNonCrossingQuadrant2
+// testShellAndHole
 template<>
 template<>
 void object::test<2> ()
 {
-    checkNonCrossing("LINESTRING (500 1000, 1000 1000, 1000 1500)",
-        "LINESTRING (300 1200, 1000 1000, 500 1500)");
+    std::vector<const NodeSection *> input{
+            sectionShell( 1,1, 5,5, 9,9 ),
+            sectionHole(  6,0, 5,5, 4,0 )
+    };
+    std::vector<const NodeSection *> expected{
+            sectionShell( 1,1, 5,5, 4,0 ),
+            sectionShell( 6,0, 5,5, 9,9 )
+   };
+    checkConversion(input, expected);
 }
 
-// testNonCrossingQuadrant4
+// testShellsAndHoles
 template<>
 template<>
 void object::test<3> ()
 {
-    checkNonCrossing("LINESTRING (500 1000, 1000 1000, 1000 1500)",
-        "LINESTRING (1000 500, 1000 1000, 1500 1000)");
+    std::vector<const NodeSection *> input{
+        sectionShell( 1,1, 5,5, 9,9 ),
+        sectionHole(  6,0, 5,5, 4,0 ),
+        
+        sectionShell( 8,8, 5,5, 1,8 ),
+        sectionHole(  4,8, 5,5, 6,8 )
+    };
+    std::vector<const NodeSection *> expected{
+        sectionShell( 1,1, 5,5, 4,0 ),
+        sectionShell( 6,0, 5,5, 9,9 ),
+        
+        sectionShell( 4,8, 5,5, 1,8 ),
+        sectionShell( 8,8, 5,5, 6,8 )
+   };
+    checkConversion(input, expected);
 }
 
-// testNonCrossingCollinear
-template<>
-template<>
-void object::test<4> ()
-{
-    checkNonCrossing("LINESTRING (3 1, 5 5, 9 9)",
-        "LINESTRING (2 1, 5 5, 9 9)");
-}
-
-// testNonCrossingBothCollinear
+// testShellAnd2Holes
 template<>
 template<>
 void object::test<5> ()
 {
-    checkNonCrossing("LINESTRING (3 1, 5 5, 9 9)",
-        "LINESTRING (3 1, 5 5, 9 9)");
+    std::vector<const NodeSection *> input{
+        sectionShell( 1,1, 5,5, 9,9 ),
+        sectionHole(  7,0, 5,5, 6,0 ),
+        sectionHole(  4,0, 5,5, 3,0 )
+    };
+    std::vector<const NodeSection *> expected{
+        sectionShell( 1,1, 5,5, 3,0 ),
+        sectionShell( 4,0, 5,5, 6,0 ),
+        sectionShell( 7,0, 5,5, 9,9 )
+    };
+    checkConversion(input, expected);
 }
 
-// testInteriorSegment
+// testHoles
 template<>
 template<>
 void object::test<6> ()
 {
-    checkInterior("LINESTRING (5 9, 5 5, 9 5)",
-        "LINESTRING (5 5, 0 0)");
+    std::vector<const NodeSection *> input{
+        sectionHole(  7,0, 5,5, 6,0 ),
+        sectionHole(  4,0, 5,5, 3,0 )
+    };
+    std::vector<const NodeSection *> expected{
+        sectionShell( 4,0, 5,5, 6,0 ),
+        sectionShell( 7,0, 5,5, 3,0 )
+    };
+    checkConversion(input, expected);
 }
 
-// testExteriorSegment
-template<>
-template<>
-void object::test<7> ()
-{
-    checkExterior("LINESTRING (5 9, 5 5, 9 5)",
-        "LINESTRING (5 5, 9 9)");
-}
-
-
 } // namespace tut
diff --git a/tests/unit/operation/relateng/RelateGeometryTest.cpp b/tests/unit/operation/relateng/RelateGeometryTest.cpp
new file mode 100644
index 000000000..b3718ec3a
--- /dev/null
+++ b/tests/unit/operation/relateng/RelateGeometryTest.cpp
@@ -0,0 +1,133 @@
+//
+// Test Suite for geos::operation::relateng::RelateGeometry 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/geom/Location.h>
+#include <geos/operation/relateng/RelateGeometry.h>
+#include <geos/operation/relateng/DimensionLocation.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_relategeometry_data {
+
+    WKTReader r;
+    // WKTWriter w;
+
+    void checkDimension(const std::string& wkt, int expectedDim, int expectedDimReal)
+    {
+        std::unique_ptr<Geometry> geom = r.read(wkt);
+        RelateGeometry rgeom(geom.get());
+        ensure_equals("checkNodeLocation 1", rgeom.getDimension(), expectedDim);
+        ensure_equals("checkNodeLocation 2", rgeom.getDimensionReal(), expectedDimReal);
+    }
+
+};
+
+typedef test_group<test_relategeometry_data> group;
+typedef group::object object;
+
+group test_relategeometry_group("geos::operation::relateng::RelateGeometry");
+
+
+// testUniquePoints
+template<>
+template<>
+void object::test<1> ()
+{
+    std::unique_ptr<Geometry> geom = r.read("MULTIPOINT ((0 0), (5 5), (5 0), (0 0))");
+    RelateGeometry rgeom(geom.get());
+    auto pts = rgeom.getUniquePoints();
+    ensure_equals("Unique pts size", pts.size(), 3u);
+}
+
+// testBoundary
+template<>
+template<>
+void object::test<2> ()
+{
+    std::unique_ptr<Geometry> geom = r.read("MULTILINESTRING ((0 0, 9 9), (9 9, 5 1))");
+    RelateGeometry rgeom(geom.get());
+    ensure("hasBoundary", rgeom.hasBoundary());
+}
+
+// testHasDimension
+template<>
+template<>
+void object::test<3> ()
+{
+    std::unique_ptr<Geometry> geom = r.read("GEOMETRYCOLLECTION (POLYGON ((1 9, 5 9, 5 5, 1 5, 1 9)), LINESTRING (1 1, 5 4), POINT (6 5))");
+    RelateGeometry rgeom(geom.get());
+    ensure("hasDimension 0", rgeom.hasDimension(0));
+    ensure("hasDimension 1", rgeom.hasDimension(1));
+    ensure("hasDimension 2", rgeom.hasDimension(2));
+}
+
+// testDimension
+template<>
+template<>
+void object::test<4> ()
+{
+    checkDimension("POINT (0 0)", 0, 0);
+}
+
+template<>
+template<>
+void object::test<5> ()
+{
+    checkDimension("LINESTRING (0 0, 0 0)", 1, 0);
+}
+
+template<>
+template<>
+void object::test<6> ()
+{
+    checkDimension("LINESTRING (0 0, 9 9)", 1, 1);
+}
+
+template<>
+template<>
+void object::test<7> ()
+{
+    checkDimension("POLYGON ((1 9, 5 9, 5 5, 1 5, 1 9))", 2, 2);
+}
+
+template<>
+template<>
+void object::test<8> ()
+{
+    checkDimension("GEOMETRYCOLLECTION (POLYGON ((1 9, 5 9, 5 5, 1 5, 1 9)), LINESTRING (1 1, 5 4), POINT (6 5))", 2, 2);
+}
+
+template<>
+template<>
+void object::test<9> ()
+{
+    checkDimension("GEOMETRYCOLLECTION (POLYGON EMPTY, LINESTRING (1 1, 5 4), POINT (6 5))", 2, 1);
+}
+
+
+
+
+
+
+
+
+} // namespace tut
diff --git a/tests/unit/operation/relateng/RelatePredicateTest.cpp b/tests/unit/operation/relateng/RelatePredicateTest.cpp
new file mode 100644
index 000000000..7408f3275
--- /dev/null
+++ b/tests/unit/operation/relateng/RelatePredicateTest.cpp
@@ -0,0 +1,138 @@
+//
+// Test Suite for geos::operation::relateng::LinearBoundary class.
+
+#include <tut/tut.hpp>
+#include <utility.h>
+
+// geos
+#include <geos/geom/Location.h>
+#include <geos/geom/Dimension.h>
+#include <geos/operation/relateng/TopologyPredicate.h>
+#include <geos/operation/relateng/RelatePredicate.h>
+
+// std
+#include <memory>
+
+
+using namespace geos::geom;
+using namespace geos::operation::relateng;
+
+
+namespace tut {
+//
+// Test Group
+//
+
+// Common data used by all tests
+struct test_relatepredicate_data {
+
+    void
+    checkPredicate(TopologyPredicate& pred, const std::string& im, bool expected)
+    {
+        applyIM(im, pred);
+        checkPred(pred, expected);
+    }
+
+    void
+    checkPredicatePartial(TopologyPredicate& pred, const std::string& im, bool expected)
+    {
+        applyIM(im, pred);
+        bool isKnown = pred.isKnown();
+        ensure("predicate value is not known", isKnown);
+        checkPred(pred, expected);
+    }
+
+    void
+    checkPred(TopologyPredicate& pred, bool expected)
+    {
+        pred.finish();
+        bool actual = pred.value();
+        ensure_equals(expected, actual);
+    }
+
+    static void
+    applyIM(const std::string& imIn, TopologyPredicate& pred)
+    {
+        std::array<Location,3> locs = {
+            Location::INTERIOR, Location::BOUNDARY, Location::EXTERIOR };
+
+        const std::string& im = cleanIM(imIn);
+        uint32_t i = 0;
+        for (Location locA : locs) {
+            for (Location locB : locs) {
+                char entry = im[i++];
+                if (entry == '0' || entry == '1' || entry == '2') {
+                    int dim = Dimension::toDimensionValue(entry);
+                    pred.updateDimension(locA, locB, dim);
+                }
+            }
+        }
+    }
+
+    static std::string
+    cleanIM(const std::string& im)
+    {
+        std::string str = im;
+        str.erase(std::remove(str.begin(), str.end(), '.'), str.end());
+        return str;
+    }
+
+    const std::string A_EXT_B_INT = "***.***.1**";
+    const std::string A_INT_B_INT = "1**.***.***";
+
+
+};
+
+typedef test_group<test_relatepredicate_data> group;
+typedef group::object object;
+
+group test_relatepredicate_group("geos::operation::relateng::RelatePredicate");
+
+//
+// Test Cases
+//
+
+// testIntersects
+template<>
+template<>
+void object::test<1> ()
+{
+    checkPredicate(*RelatePredicate::intersects(), A_INT_B_INT, true);
+}
+
+// testDisjoint
+template<>
+template<>
+void object::test<2> ()
+{
+    checkPredicate(*RelatePredicate::intersects(), A_EXT_B_INT, false);
+    checkPredicate(*RelatePredicate::disjoint(), A_EXT_B_INT, true);
+}
+
+// testCovers
+template<>
+template<>
+void object::test<3> ()
+{
+    checkPredicate(*RelatePredicate::covers(), A_INT_B_INT, true);
+    checkPredicate(*RelatePredicate::covers(), A_EXT_B_INT, false);
+}
+
+// testCoversFast
+template<>
+template<>
+void object::test<4> ()
+{
+    checkPredicatePartial(*RelatePredicate::covers(), A_EXT_B_INT, false);
+}
+
+// testMatch
+template<>
+template<>
+void object::test<5> ()
+{
+    checkPredicate(*RelatePredicate::matches("1***T*0**"), "1**.*2*.0**", true);
+}
+
+
+} // namespace tut

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

Summary of changes:
 NEWS.md                                            |  19 +-
 benchmarks/capi/CMakeLists.txt                     |   2 +-
 capi/geos_c.cpp                                    |   2 +-
 capi/geos_c.h.in                                   |  10 +-
 include/geos/algorithm/ConvexHull.h                |   4 +-
 include/geos/algorithm/HCoordinate.h               |   2 +-
 include/geos/algorithm/LineIntersector.h           |  16 +-
 include/geos/constants.h                           |   1 +
 include/geos/coverage/CoverageRingEdges.h          |   1 +
 include/geos/geom/CoordinateList.h                 |   2 +-
 include/geos/geom/CoordinateSequences.h            |   2 +-
 include/geos/geom/Geometry.h                       |   1 +
 include/geos/geom/GeometryCollection.h             |   8 +
 include/geos/geom/GeometryFactory.h                |   4 +-
 include/geos/geom/IntersectionMatrix.h             |   5 +-
 include/geos/geom/PrecisionModel.h                 |   2 +-
 include/geos/geom/util/GeometryEditor.h            |   2 +-
 include/geos/geom/util/GeometryTransformer.h       |   2 +-
 include/geos/geomgraph/DirectedEdge.h              |   4 +-
 include/geos/geomgraph/DirectedEdgeStar.h          |   2 +-
 include/geos/index/bintree/Bintree.h               |   2 +-
 include/geos/index/quadtree/Quadtree.h             |   2 +-
 include/geos/io/GeoJSONReader.h                    |   4 +-
 include/geos/io/WKBReader.h                        |   2 +-
 include/geos/linearref/LengthIndexedLine.h         |   2 +-
 include/geos/namespaces.h                          |   2 +-
 include/geos/noding/Octant.h                       |   2 +-
 include/geos/noding/SegmentString.h                |  36 ++
 include/geos/operation/buffer/BufferOp.h           |   1 +
 include/geos/operation/buffer/OffsetCurveSection.h |   1 +
 .../geos/operation/overlayng/EdgeNodingBuilder.h   |   2 +-
 include/geos/operation/polygonize/BuildArea.h      |   2 +-
 .../geos/operation/polygonize/PolygonizeGraph.h    |   2 +-
 include/geos/operation/relate/EdgeEndBundle.h      |   2 +-
 .../geos/operation/relateng/AdjacentEdgeLocator.h  |   2 +-
 include/geos/operation/relateng/BasicPredicate.h   | 105 ++++
 include/geos/operation/relateng/IMPatternMatcher.h |  87 +++
 include/geos/operation/relateng/IMPredicate.h      | 131 +++++
 .../operation/relateng/IntersectionMatrixPattern.h |  65 +++
 include/geos/operation/relateng/RelateGeometry.h   | 161 +++++-
 .../operation/relateng/RelateMatrixPredicate.h     |  85 +++
 include/geos/operation/relateng/RelatePredicate.h  | 639 +++++++++++++++++++++
 .../geos/operation/relateng/RelateSegmentString.h  | 170 ++++++
 .../geos/operation/relateng/TopologyPredicate.h    | 199 +++++++
 include/geos/operation/sharedpaths/SharedPathsOp.h |   8 +-
 .../geos/operation/union/CascadedPolygonUnion.h    |   2 +-
 .../operation/valid/IndexedNestedPolygonTester.h   |   2 +-
 include/geos/planargraph/DirectedEdgeStar.h        |   2 +-
 include/geos/planargraph/Subgraph.h                |   2 +-
 include/geos/precision/CommonBitsOp.h              |   2 +-
 include/geos/precision/GeometryPrecisionReducer.h  |   2 +-
 include/geos/simplify/TaggedLineSegment.h          |   2 +-
 .../triangulate/quadedge/QuadEdgeSubdivision.h     |   2 +-
 include/geos/vend/json.hpp                         |  14 +-
 src/algorithm/CGAlgorithmsDD.cpp                   |   2 +-
 src/algorithm/RobustDeterminant.cpp                |   2 +-
 src/geom/CoordinateSequence.cpp                    |   2 +-
 src/geom/GeometryCollection.cpp                    |  19 +
 src/operation/buffer/BufferBuilder.cpp             |   2 +-
 src/operation/buffer/BufferCurveSetBuilder.cpp     |   2 +-
 src/operation/buffer/BufferSubgraph.cpp            |   2 +-
 src/operation/buffer/OffsetCurve.cpp               |   6 +-
 src/operation/buffer/SubgraphDepthLocater.cpp      |   2 +-
 src/operation/linemerge/LineSequencer.cpp          |   2 +-
 src/operation/overlayng/OverlayGraph.cpp           |   2 +-
 src/operation/polygonize/BuildArea.cpp             |   2 +-
 src/operation/relateng/BasicPredicate.cpp          | 137 +++++
 src/operation/relateng/IMPatternMatcher.cpp        | 149 +++++
 src/operation/relateng/IMPredicate.cpp             | 155 +++++
 src/operation/relateng/LinearBoundary.cpp          |   3 -
 src/operation/relateng/NodeSection.cpp             |   8 +-
 src/operation/relateng/RelateEdge.cpp              |   1 -
 src/operation/relateng/RelateGeometry.cpp          | 464 +++++++++++++++
 src/operation/relateng/RelatePredicate.cpp         | 109 ++++
 src/operation/relateng/RelateSegmentString.cpp     | 224 ++++++++
 src/operation/sharedpaths/SharedPathsOp.cpp        |   2 +-
 src/simplify/LinkedLine.cpp                        |   2 +-
 tests/README.md                                    |   2 +-
 .../distance/DiscreteFrechetDistanceTest.cpp       |   2 +-
 tests/unit/capi/GEOSGeom_createRectangleTest.cpp   |   2 +-
 tests/unit/capi/GEOSIntersectionTest.cpp           |  16 +
 tests/unit/geom/CoordinateSequenceTest.cpp         |   6 +-
 .../unit/operation/relateng/RelateGeometryTest.cpp | 133 +++++
 .../operation/relateng/RelatePredicateTest.cpp     | 138 +++++
 tests/unit/tut/tut.hpp                             |   4 +-
 tests/unit/tut/tut_exception.hpp                   |  10 +-
 tests/unit/tut/tut_result.hpp                      |   6 +-
 tests/xmltester/XMLTester.cpp                      |   6 +-
 tests/xmltester/tests/misc/InvalidRelates.xml      |   2 +-
 tests/xmltester/tests/validate/TestRelatePA.xml    |   2 +-
 tests/xmltester/tinyxml2/tinyxml2.cpp              |   2 +-
 util/geosop/GeometryOp.cpp                         |   2 +-
 web/content/project/psc.md                         |   6 +-
 web/content/project/rfcs/rfc01.md                  |   2 +-
 web/content/project/rfcs/rfc06.md                  |   2 +-
 web/content/specifications/wkb.md                  |   2 +-
 web/content/usage/faq/index.md                     |   2 +-
 web/content/usage/install.md                       |   2 +-
 98 files changed, 3362 insertions(+), 121 deletions(-)
 create mode 100644 include/geos/operation/relateng/BasicPredicate.h
 create mode 100644 include/geos/operation/relateng/IMPatternMatcher.h
 create mode 100644 include/geos/operation/relateng/IMPredicate.h
 create mode 100644 include/geos/operation/relateng/IntersectionMatrixPattern.h
 create mode 100644 include/geos/operation/relateng/RelateMatrixPredicate.h
 create mode 100644 include/geos/operation/relateng/RelatePredicate.h
 create mode 100644 include/geos/operation/relateng/RelateSegmentString.h
 create mode 100644 include/geos/operation/relateng/TopologyPredicate.h
 create mode 100644 src/operation/relateng/BasicPredicate.cpp
 create mode 100644 src/operation/relateng/IMPatternMatcher.cpp
 create mode 100644 src/operation/relateng/IMPredicate.cpp
 create mode 100644 src/operation/relateng/RelateGeometry.cpp
 create mode 100644 src/operation/relateng/RelatePredicate.cpp
 create mode 100644 src/operation/relateng/RelateSegmentString.cpp
 create mode 100644 tests/unit/operation/relateng/RelateGeometryTest.cpp
 create mode 100644 tests/unit/operation/relateng/RelatePredicateTest.cpp


hooks/post-receive
-- 
GEOS