Fix RelateOp for empty geometry and closed linear geometry

Ports the RelateOp class from JTS and then applies the fixes
in locationtech/jts#671

Fixes https://trac.osgeo.org/geos/ticket/1096

include/geos/operation/BoundaryOp.h

@@ -0,0 +1,120 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2022 ISciences LLC
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************
+ *
+ * Last port: operation/BoundaryOp.java fd5aebb
+ *
+ **********************************************************************/
+
+#pragma once
+
+#include <geos/algorithm/BoundaryNodeRule.h>
+#include <geos/geom/Geometry.h>
+#include <map>
+
+namespace geos {
+namespace geom {
+class LineString;
+class MultiLineString;
+}
+}
+
+namespace geos {
+namespace operation {
+
+/**
+ * Computes the boundary of a Geometry.
+ * Allows specifying the BoundaryNodeRule to be used.
+ * This operation will always return a Geometry of the appropriate
+ * dimension for the boundary (even if the input geometry is empty).
+ * The boundary of zero-dimensional geometries (Points) is
+ * always the empty GeometryCollection.
+ *
+ * @author Martin Davis
+ * @version 1.7
+ */
+class GEOS_DLL BoundaryOp {
+
+public:
+  /**
+   * Creates a new instance for the given geometry.
+   *
+   * @param geom the input geometry
+   */
+  BoundaryOp(const geom::Geometry& geom);
+
+  /**
+   * Creates a new instance for the given geometry.
+   *
+   * @param geom the input geometry
+   * @param bnRule the Boundary Node Rule to use
+   */
+  BoundaryOp(const geom::Geometry& geom, const algorithm::BoundaryNodeRule& bnRule);
+
+ /**
+   * Computes a geometry representing the boundary of a geometry.
+   *
+   * @param g the input geometry
+   * @return the computed boundary
+   */
+  static std::unique_ptr<geom::Geometry> getBoundary(const geom::Geometry& g);
+
+  /**
+   * Computes a geometry representing the boundary of a geometry,
+   * using an explicit BoundaryNodeRule.
+   *
+   * @param g the input geometry
+   * @param bnRule the Boundary Node Rule to use
+   * @return the computed boundary
+   */
+  static std::unique_ptr<geom::Geometry> getBoundary(const geom::Geometry& g, const algorithm::BoundaryNodeRule& bnRule);
+
+  /**
+   * Tests if a geometry has a boundary (it is non-empty).
+   * The semantics are:
+   * <ul>
+   * <li>Empty geometries do not have boundaries.
+   * <li>Points do not have boundaries.
+   * <li>For linear geometries the existence of the boundary
+   * is determined by the BoundaryNodeRule.
+   * <li>Non-empty polygons always have a boundary.
+   * </ul>
+   *
+   * @param geom the geometry providing the boundary
+   * @param boundaryNodeRule  the Boundary Node Rule to use
+   * @return true if the boundary exists
+   */
+  static bool hasBoundary(const geom::Geometry& geom, const algorithm::BoundaryNodeRule& boundaryNodeRule);
+
+  /**
+   * Gets the computed boundary.
+   *
+   * @return the boundary geometry
+   */
+  std::unique_ptr<geom::Geometry> getBoundary();
+
+private:
+  const geom::Geometry& m_geom;
+  const geom::GeometryFactory& m_geomFact;
+  const algorithm::BoundaryNodeRule& m_bnRule;
+
+  std::unique_ptr<geom::Geometry> boundaryMultiLineString(const geom::MultiLineString& mLine);
+
+  std::vector<geom::Coordinate> computeBoundaryCoordinates(const geom::MultiLineString& mLine);
+
+  std::unique_ptr<geom::Geometry> boundaryLineString(const geom::LineString& line);
+};
+
+}
+}
+

include/geos/operation/relate/RelateComputer.h

@@ -38,6 +38,9 @@
 
 // Forward declarations
 namespace geos {
+namespace algorithm {
+class BoundaryNodeRule;
+}
 namespace geom {
 class Geometry;
 }
@@ -110,7 +113,8 @@ class GEOS_DLL RelateComputer {
      * If the Geometries are disjoint, we need to enter their dimension and
      * boundary dimension in the Ext rows in the IM
      */
-    void computeDisjointIM(geom::IntersectionMatrix* imX);
+    void computeDisjointIM(geom::IntersectionMatrix* imX,
+                           const algorithm::BoundaryNodeRule& boundaryNodeRule);
 
     void labelNodeEdges();
 
@@ -119,6 +123,21 @@ class GEOS_DLL RelateComputer {
      */
     void updateIM(geom::IntersectionMatrix& imX);
 
+    /**
+     * Compute the IM entry for the intersection of the boundary
+     * of a geometry with the Exterior.
+     * This is the nominal dimension of the boundary
+     * unless the boundary is empty, in which case it is {@link Dimension#FALSE}.
+     * For linear geometries the Boundary Node Rule determines
+     * whether the boundary is empty.
+     *
+     * @param geom the geometry providing the boundary
+     * @param boundaryNodeRule  the Boundary Node Rule to use
+     * @return the IM dimension entry
+     */
+    static int getBoundaryDim(const geom::Geometry& geom,
+                              const algorithm::BoundaryNodeRule& boundaryNodeRule);
+
     /**
      * Processes isolated edges by computing their labelling and adding them
      * to the isolated edges list.

src/geom/LineString.cpp

@@ -35,6 +35,7 @@
 #include <geos/geom/Point.h>
 #include <geos/geom/MultiPoint.h> // for getBoundary
 #include <geos/geom/Envelope.h>
+#include <geos/operation/BoundaryOp.h>
 
 #include <algorithm>
 #include <typeinfo>
@@ -233,19 +234,8 @@ LineString::getGeometryType() const
 std::unique_ptr<Geometry>
 LineString::getBoundary() const
 {
-    if(isEmpty()) {
-        return std::unique_ptr<Geometry>(getFactory()->createMultiPoint());
-    }
-
-    // using the default OGC_SFS MOD2 rule, the boundary of a
-    // closed LineString is empty
-    if(isClosed()) {
-        return std::unique_ptr<Geometry>(getFactory()->createMultiPoint());
-    }
-    std::vector<std::unique_ptr<Point>> pts(2);
-    pts[0] = getStartPoint();
-    pts[1] = getEndPoint();
-    return getFactory()->createMultiPoint(std::move(pts));
+    operation::BoundaryOp bop(*this);
+    return bop.getBoundary();
 }
 
 bool

src/geom/MultiLineString.cpp

@@ -21,6 +21,7 @@
 #include <geos/geomgraph/GeometryGraph.h>
 #include <geos/geom/GeometryFactory.h>
 #include <geos/geom/MultiLineString.h>
+#include <geos/operation/BoundaryOp.h>
 
 #include <vector>
 #include <cassert>
@@ -89,13 +90,8 @@ MultiLineString::isClosed() const
 std::unique_ptr<Geometry>
 MultiLineString::getBoundary() const
 {
-    if(isEmpty()) {
-        return getFactory()->createGeometryCollection();
-    }
-
-    GeometryGraph gg(0, this);
-    CoordinateSequence* pts = gg.getBoundaryPoints();
-    return std::unique_ptr<Geometry>(getFactory()->createMultiPoint(*pts));
+    operation::BoundaryOp bop(*this);
+    return bop.getBoundary();
 }
 
 GeometryTypeId

src/operation/BoundaryOp.cpp

@@ -0,0 +1,173 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2022 ISciences LLC
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************
+ *
+ * Last port: operation/BoundaryOp.java fd5aebb
+ *
+ **********************************************************************/
+
+#include <geos/operation/BoundaryOp.h>
+#include <geos/algorithm/BoundaryNodeRule.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/GeometryFactory.h>
+#include <geos/geom/LineString.h>
+#include <geos/geom/MultiLineString.h>
+#include <map>
+
+using geos::geom::Coordinate;
+using geos::geom::Dimension;
+using geos::geom::Geometry;
+using geos::geom::LineString;
+using geos::geom::MultiLineString;
+using geos::geom::Point;
+using geos::algorithm::BoundaryNodeRule;
+
+namespace geos {
+namespace operation {
+
+BoundaryOp::BoundaryOp(const Geometry& geom) :
+    m_geom(geom),
+    m_geomFact(*geom.getFactory()),
+    m_bnRule(BoundaryNodeRule::getBoundaryRuleMod2())
+{}
+
+BoundaryOp::BoundaryOp(const geom::Geometry& geom, const algorithm::BoundaryNodeRule& bnRule) :
+    m_geom(geom),
+    m_geomFact(*geom.getFactory()),
+    m_bnRule(bnRule)
+{}
+
+std::unique_ptr<geom::Geometry>
+BoundaryOp::getBoundary()
+{
+    if (auto ls = dynamic_cast<const LineString*>(&m_geom)) {
+        return boundaryLineString(*ls);
+    }
+
+    if (auto mls = dynamic_cast<const MultiLineString*>(&m_geom)) {
+        return boundaryMultiLineString(*mls);
+    }
+
+    return m_geom.getBoundary();
+}
+
+std::unique_ptr<geom::Geometry>
+BoundaryOp::getBoundary(const geom::Geometry& g)
+{
+    BoundaryOp bop(g);
+    return bop.getBoundary();
+}
+
+std::unique_ptr<geom::Geometry>
+BoundaryOp::getBoundary(const geom::Geometry& g, const algorithm::BoundaryNodeRule& bnRule)
+{
+    BoundaryOp bop(g, bnRule);
+    return bop.getBoundary();
+}
+
+bool
+BoundaryOp::hasBoundary(const geom::Geometry& geom, const algorithm::BoundaryNodeRule& boundaryNodeRule)
+{
+    // Note that this does not handle geometry collections with a non-empty linear element
+    if (geom.isEmpty()) {
+        return false;
+    }
+
+    switch (geom.getDimension()) {
+        case Dimension::P: return false;
+        /**
+         * Linear geometries might have an empty boundary due to boundary node rule.
+         */
+        case Dimension::L:
+            {
+
+            auto boundary = getBoundary(geom, boundaryNodeRule);
+            return !boundary->isEmpty();
+            }
+        default:
+            return true;
+        }
+}
+
+std::unique_ptr<Geometry>
+BoundaryOp::boundaryLineString(const geom::LineString& line)
+{
+    if (m_geom.isEmpty()) {
+        return m_geomFact.createMultiPoint();
+    }
+
+    if (line.isClosed()) {
+        // check whether endpoints of valence 2 are on the boundary or not
+        bool closedEndpointOnBoundary = m_bnRule.isInBoundary(2);
+        if (closedEndpointOnBoundary) {
+            return line.getStartPoint();
+        }
+        else {
+            return m_geomFact.createMultiPoint();
+        }
+    }
+
+    std::vector<std::unique_ptr<Point>> pts(2);
+    pts[0] = line.getStartPoint();
+    pts[1] = line.getEndPoint();
+
+    return m_geomFact.createMultiPoint(std::move(pts));
+}
+
+std::unique_ptr<Geometry>
+BoundaryOp::boundaryMultiLineString(const geom::MultiLineString& mLine)
+{
+    if (m_geom.isEmpty()) {
+        return m_geomFact.createMultiPoint();
+    }
+
+    auto bdyPts = computeBoundaryCoordinates(mLine);
+
+    // return Point or MultiPoint
+    if (bdyPts.size() == 1) {
+        return std::unique_ptr<Geometry>(m_geomFact.createPoint(bdyPts[0]));
+    }
+    // this handles 0 points case as well
+    return m_geomFact.createMultiPoint(std::move(bdyPts));
+}
+
+std::vector<geom::Coordinate>
+BoundaryOp::computeBoundaryCoordinates(const geom::MultiLineString& mLine)
+{
+    std::vector<Coordinate> bdyPts;
+    std::map<Coordinate, int> endpointMap;
+
+    for (std::size_t i = 0; i < mLine.getNumGeometries(); i++) {
+      const LineString* line = mLine.getGeometryN(i);
+
+      if (line->getNumPoints() == 0) {
+        continue;
+      }
+
+      endpointMap[line->getCoordinateN(0)]++;
+      endpointMap[line->getCoordinateN(line->getNumPoints() - 1)]++;
+    }
+
+    for (const auto& entry: endpointMap) {
+        auto valence = entry.second;
+        if (m_bnRule.isInBoundary(valence)) {
+            bdyPts.push_back(entry.first);
+        }
+    }
+
+    return bdyPts;
+}
+
+
+}
+}