Add GEOSOrientPolygons (libgeos#818)

* Add GEOSOrientPolygons

Fixes libgeos#779

Co-authored-by: Joris Van den Bossche <[email protected]>

NEWS.md

@@ -15,6 +15,7 @@ xxxx-xx-xx
   - CAPI: GEOSDisjointSubsetUnion (GH-692, Dan Baston)
   - CAPI: GEOSLineSubstring (GH-706, Dan Baston)
   - CAPI: GEOSEqualsIdentical (GH-810, Dan Baston)
+  - CAPI: GEOSOrientPolygons (GH-818, Dan Baston)
   - CAPI: GEOSSTRtree_build (GH-835, Dan Baston)
   - CAPI: GEOSConcaveHullByLength (GH-849, Martin Davis)
   - CAPI: GEOSGeomGetM (GH-864, Mike Taves)

capi/geos_c.cpp

@@ -681,6 +681,12 @@ extern "C" {
         return GEOSNormalize_r(handle, g);
     }
 
+    int
+    GEOSOrientPolygons(Geometry* g, int exteriorCW)
+    {
+        return GEOSOrientPolygons_r(handle, g, exteriorCW);
+    }
+
     int
     GEOSGetNumInteriorRings(const Geometry* g)
     {

capi/geos_c.h.in

@@ -1535,6 +1535,12 @@ extern int GEOS_DLL GEOSNormalize_r(
     GEOSContextHandle_t handle,
     GEOSGeometry* g);
 
+/** \see GEOSOrientPolygons */
+extern int GEOS_DLL GEOSOrientPolygons_r(
+    GEOSContextHandle_t handle,
+    GEOSGeometry* g,
+    int exteriorCW);
+
 /**
 * Controls the behavior of GEOSGeom_setPrecision()
 * when altering the precision of a geometry.
@@ -2916,6 +2922,18 @@ extern void GEOS_DLL GEOSGeom_setUserData(GEOSGeometry* g, void* userData);
 */
 extern int GEOS_DLL GEOSNormalize(GEOSGeometry* g);
 
+/**
+* Enforce a ring orientation on all polygonal elements in the input geometry.
+* Non-polygonal geometries will not be modified.
+*
+* \param g Input geometry
+* \param exteriorCW if 1, exterior rings will be clockwise and interior rings
+*                         will be counter-clockwise
+* \return 0 on success or -1 on exception
+*/
+extern int GEOS_DLL GEOSOrientPolygons(GEOSGeometry* g,
+                                       int exteriorCW);
+
 ///@}
 
 /* ========== Validity checking ============================================================ */

capi/geos_ts_c.cpp

@@ -1671,6 +1671,32 @@ extern "C" {
         });
     }
 
+    int
+    GEOSOrientPolygons_r(GEOSContextHandle_t extHandle, Geometry* g, int exteriorCW)
+    {
+        return execute(extHandle, -1, [&]() {
+            class OrientPolygons : public geos::geom::GeometryComponentFilter {
+            public:
+                OrientPolygons(bool isExteriorCW) : exteriorCW(isExteriorCW) {}
+
+                void filter_rw(Geometry* g) override {
+                    if (g->getGeometryTypeId() == geos::geom::GeometryTypeId::GEOS_POLYGON) {
+                        auto p = geos::detail::down_cast<Polygon*>(g);
+                        p->orientRings(exteriorCW);
+                    }
+                }
+
+            private:
+                bool exteriorCW;
+            };
+
+            OrientPolygons op(exteriorCW);
+            g->apply_rw(&op);
+
+            return 0;
+        });
+    }
+
     int
     GEOSGetNumInteriorRings_r(GEOSContextHandle_t extHandle, const Geometry* g1)
     {

include/geos/geom/LinearRing.h

@@ -103,6 +103,8 @@ class GEOS_DLL LinearRing : public LineString {
 
     std::unique_ptr<LinearRing> reverse() const { return std::unique_ptr<LinearRing>(reverseImpl()); }
 
+    void orient(bool isCW);
+
 protected:
 
     int

include/geos/geom/Polygon.h

@@ -152,6 +152,16 @@ class GEOS_DLL Polygon: public Geometry {
 
     void normalize() override;
 
+    /**
+     * \brief
+     * Apply a ring ordering convention to this polygon, with
+     * interior rings having an opposite orientation to the
+     * specified exterior orientation.
+     *
+     * \param exteriorCW should exterior ring be clockwise?
+     */
+    void orientRings(bool exteriorCW);
+
     std::unique_ptr<Polygon> reverse() const { return std::unique_ptr<Polygon>(reverseImpl()); }
 
     const CoordinateXY* getCoordinate() const override;

src/geom/LinearRing.cpp

@@ -17,6 +17,8 @@
  *
  **********************************************************************/
 
+#include <geos/algorithm/Orientation.h>
+
 #include <geos/geom/LinearRing.h>
 #include <geos/geom/Dimension.h>
 #include <geos/geom/CoordinateSequence.h>
@@ -98,6 +100,19 @@ LinearRing::getGeometryTypeId() const
     return GEOS_LINEARRING;
 }
 
+void
+LinearRing::orient(bool isCW)
+{
+    if (isEmpty()) {
+        return;
+    }
+
+    if (algorithm::Orientation::isCCW(points.get()) == isCW) {
+        points->reverse();
+    }
+
+}
+
 LinearRing*
 LinearRing::reverseImpl() const
 {

src/geom/Polygon.cpp

@@ -363,6 +363,15 @@ Polygon::normalize()
     });
 }
 
+void
+Polygon::orientRings(bool exteriorCW)
+{
+    shell->orient(exteriorCW);
+    for (auto& hole : holes) {
+        hole->orient(!exteriorCW);
+    }
+}
+
 int
 Polygon::compareToSameClass(const Geometry* g) const
 {

tests/unit/capi/GEOSOrientPolygonsTest.cpp

@@ -0,0 +1,94 @@
+//
+// Test Suite for C-API GEOSOrientationPolygons
+
+#include <tut/tut.hpp>
+// geos
+#include <geos_c.h>
+// std
+#include <string>
+#include <cstdarg>
+#include <cstdio>
+#include <cstdlib>
+#include <memory>
+
+#include "capi_test_utils.h"
+
+namespace tut {
+//
+// Test Group
+//
+
+// Common data used in test cases.
+struct test_capigeosorientpolygons_data : public capitest::utility {};
+
+typedef test_group<test_capigeosorientpolygons_data> group;
+typedef group::object object;
+
+group test_capigeosorientpolygons_group("capi::GEOSOrientPolygons");
+
+// empty polygon is passed through
+template<>
+template<>
+void object::test<1>()
+{
+    geom1_ = GEOSGeomFromWKT("POLYGON EMPTY");
+    ensure_equals(GEOSOrientPolygons(geom1_, 1), 0);
+
+    ensure_equals(toWKT(geom1_), "POLYGON EMPTY");
+}
+
+// hole orientation is opposite to shell
+template<>
+template<>
+void object::test<2>()
+{
+    geom1_ = GEOSGeomFromWKT("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0), (1 1, 2 1, 2 2, 1 2, 1 1))");
+
+    ensure_equals(GEOSOrientPolygons(geom1_, 0), 0);
+    ensure_equals(toWKT(geom1_), "POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0), (1 1, 1 2, 2 2, 2 1, 1 1))");
+
+    ensure_equals(GEOSOrientPolygons(geom1_, 1), 0);
+    ensure_equals(toWKT(geom1_), "POLYGON ((0 0, 0 10, 10 10, 10 0, 0 0), (1 1, 2 1, 2 2, 1 2, 1 1))");
+}
+
+// all polygons in collection are processed
+template<>
+template<>
+void object::test<3>()
+{
+    geom1_ = GEOSGeomFromWKT("MULTIPOLYGON (((0 0, 10 0, 10 10, 0 10, 0 0), (1 1, 2 1, 2 2, 1 2, 1 1)), ((100 100, 200 100, 200 200, 100 100)))");
+
+    ensure_equals(GEOSOrientPolygons(geom1_, 0), 0);
+    ensure_equals(toWKT(geom1_), "MULTIPOLYGON (((0 0, 10 0, 10 10, 0 10, 0 0), (1 1, 1 2, 2 2, 2 1, 1 1)), ((100 100, 200 100, 200 200, 100 100)))");
+
+    ensure_equals(GEOSOrientPolygons(geom1_, 1), 0);
+    ensure_equals(toWKT(geom1_), "MULTIPOLYGON (((0 0, 0 10, 10 10, 10 0, 0 0), (1 1, 2 1, 2 2, 1 2, 1 1)), ((100 100, 200 200, 200 100, 100 100)))");
+}
+
+// polygons in collection are oriented, closed linestring unchanged
+template<>
+template<>
+void object::test<4>()
+{
+    geom1_ = GEOSGeomFromWKT("GEOMETRYCOLLECTION (POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0), (1 1, 2 1, 2 2, 1 2, 1 1)), LINESTRING (100 100, 200 100, 200 200, 100 100))");
+
+    ensure_equals(GEOSOrientPolygons(geom1_, 1), 0);
+    ensure_equals(toWKT(geom1_), "GEOMETRYCOLLECTION (POLYGON ((0 0, 0 10, 10 10, 10 0, 0 0), (1 1, 2 1, 2 2, 1 2, 1 1)), LINESTRING (100 100, 200 100, 200 200, 100 100))");
+}
+
+// nested collection handled correctly
+template<>
+template<>
+void object::test<5>()
+{
+    geom1_ = GEOSGeomFromWKT("GEOMETRYCOLLECTION (GEOMETRYCOLLECTION (MULTIPOLYGON (((0 0, 10 0, 10 10, 0 10, 0 0)))))");
+
+    ensure_equals(GEOSOrientPolygons(geom1_, 0), 0);
+    ensure_equals(toWKT(geom1_), "GEOMETRYCOLLECTION (GEOMETRYCOLLECTION (MULTIPOLYGON (((0 0, 10 0, 10 10, 0 10, 0 0)))))");
+
+    ensure_equals(GEOSOrientPolygons(geom1_, 1), 0);
+    ensure_equals(toWKT(geom1_), "GEOMETRYCOLLECTION (GEOMETRYCOLLECTION (MULTIPOLYGON (((0 0, 0 10, 10 10, 10 0, 0 0)))))");
+}
+
+} // namespace tut
+