From a13c2aab335ac237b89e0748cf74e0c15c8498cc Mon Sep 17 00:00:00 2001
From: vadimcn <vadimcn@gmail.com>
Date: Wed, 22 May 2024 17:12:28 -0700
Subject: [PATCH] GEOSClipByRect: preserve Z coordinate (#1095)

For [multi]linestrings/linerings, Z coordinate is interpolated
from that of the nearest existing vertices.

For [multi]polygons, Z coordinate of points created on the boundary
is interpolated as above.  For interior points, however, it is populated
using the ElevationModel.
---
 .../intersection/RectangleIntersection.cpp    | 54 +++++++----
 tests/unit/capi/GEOSClipByRectTest.cpp        | 91 ++++++++++++++++++-
 2 files changed, 121 insertions(+), 24 deletions(-)

diff --git a/src/operation/intersection/RectangleIntersection.cpp b/src/operation/intersection/RectangleIntersection.cpp
index e9ddc31428..e3f7750aff 100644
--- a/src/operation/intersection/RectangleIntersection.cpp
+++ b/src/operation/intersection/RectangleIntersection.cpp
@@ -17,6 +17,7 @@
 #include <geos/operation/intersection/RectangleIntersection.h>
 #include <geos/operation/intersection/Rectangle.h>
 #include <geos/operation/intersection/RectangleIntersectionBuilder.h>
+#include <geos/operation/overlayng/ElevationModel.h>
 #include <geos/operation/predicate/RectangleIntersects.h>
 #include <geos/geom/GeometryFactory.h>
 #include <geos/geom/CoordinateSequence.h>
@@ -34,6 +35,7 @@
 
 using geos::operation::intersection::Rectangle;
 using geos::operation::intersection::RectangleIntersectionBuilder;
+using geos::operation::overlayng::ElevationModel;
 using namespace geos::geom;
 using namespace geos::algorithm;
 namespace geos {
@@ -62,15 +64,18 @@ different(double x1, double y1, double x2, double y2)
 
 inline
 void
-clip_one_edge(double& x1, double& y1, double x2, double y2, double limit)
+clip_one_edge(double& x1, double& y1, double& z1, double x2, double y2, double z2, double limit)
 {
     if(x2 == limit) {
         y1 = y2;
         x1 = x2;
+        z1 = z2;
     }
 
     if(x1 != x2) {
-        y1 += (y2 - y1) * (limit - x1) / (x2 - x1);
+        double fraction = (limit - x1) / (x2 - x1);
+        y1 += (y2 - y1) * fraction;
+        z1 += (z2 - z1) * fraction;
         x1 = limit;
     }
 }
@@ -86,22 +91,22 @@ clip_one_edge(double& x1, double& y1, double x2, double y2, double limit)
  */
 
 void
-clip_to_edges(double& x1, double& y1,
-              double x2, double y2,
+clip_to_edges(double& x1, double& y1, double& z1,
+              double x2, double y2, double z2,
               const Rectangle& rect)
 {
     if(x1 < rect.xmin()) {
-        clip_one_edge(x1, y1, x2, y2, rect.xmin());
+        clip_one_edge(x1, y1, z1, x2, y2, z2, rect.xmin());
     }
     else if(x1 > rect.xmax()) {
-        clip_one_edge(x1, y1, x2, y2, rect.xmax());
+        clip_one_edge(x1, y1, z1, x2, y2, z2, rect.xmax());
     }
 
     if(y1 < rect.ymin()) {
-        clip_one_edge(y1, x1, y2, x2, rect.ymin());
+        clip_one_edge(y1, x1, z1, y2, x2, z2, rect.ymin());
     }
     else if(y1 > rect.ymax()) {
-        clip_one_edge(y1, x1, y2, x2, rect.ymax());
+        clip_one_edge(y1, x1, z1, y2, x2, z2, rect.ymax());
     }
 }
 
@@ -153,6 +158,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
 
     double x0 = 0;
     double y0 = 0;
+    double z0 = 0;
     bool add_start = false;
 
     // Start iterating
@@ -164,6 +170,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
 
         double x = cs[i].x;
         double y = cs[i].y;
+        double z = cs[i].z;
         Rectangle::Position pos = rect.position(x, y);
 
         if(pos == Rectangle::Outside) {
@@ -199,12 +206,14 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
             // Establish new position
             x = cs[i].x;
             y = cs[i].y;
+            z = cs[i].z;
             pos = rect.position(x, y);
 
             // Handle all possible cases
             x0 = cs[i - 1].x;
             y0 = cs[i - 1].y;
-            clip_to_edges(x0, y0, x, y, rect);
+            z0 = cs[i - 1].z;
+            clip_to_edges(x0, y0, z0, x, y, z, rect);
 
             if(pos == Rectangle::Inside) {
                 add_start = true;	// x0,y0 must have clipped the rectangle
@@ -218,7 +227,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
                 // which will then enter the Outside section.
 
                 // Clip the other end too
-                clip_to_edges(x, y, x0, y0, rect);
+                clip_to_edges(x, y, z, x0, y0, z0, rect);
 
                 Rectangle::Position prev_pos = rect.position(x0, y0);
                 pos = rect.position(x, y);
@@ -229,8 +238,8 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
                         !Rectangle::onSameEdge(prev_pos, pos)	// discard if travels along edge
                   ) {
                     auto coords = detail::make_unique<geom::CoordinateSequence>(2u);
-                    coords->setAt(Coordinate(x0, y0), 0);
-                    coords->setAt(Coordinate(x, y), 1);
+                    coords->setAt(Coordinate(x0, y0, z0), 0);
+                    coords->setAt(Coordinate(x, y, z), 1);
                     auto line = _gf->createLineString(std::move(coords));
                     parts.add(line.release());
                 }
@@ -267,6 +276,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
             while(!go_outside && ++i < n) {
                 x = cs[i].x;
                 y = cs[i].y;
+                z = cs[i].z;
 
                 Rectangle::Position prev_pos = pos;
                 pos = rect.position(x, y);
@@ -278,7 +288,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
                     go_outside = true;
 
                     // Clip the outside point to edges
-                    clip_to_edges(x, y, cs[i - 1].x, cs[i - 1].y, rect);
+                    clip_to_edges(x, y, z, cs[i - 1].x, cs[i - 1].y, cs[i - 1].z, rect);
                     pos = rect.position(x, y);
 
                     // Does the line exit through the inside of the box?
@@ -291,7 +301,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
                     if(start_index < i - 1 || add_start || through_box) {
                         auto coords = detail::make_unique<CoordinateSequence>();
                         if(add_start) {
-                            coords->add(Coordinate(x0, y0));
+                            coords->add(Coordinate(x0, y0, z0));
                             add_start = false;
                         }
                         //line->addSubLineString(&g, start_index, i-1);
@@ -299,7 +309,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
                                     cs.begin() + static_cast<long>(i));
 
                         if(through_box) {
-                            coords->add(Coordinate(x, y));
+                            coords->add(Coordinate(x, y, z));
                         }
 
                         auto line = _gf->createLineString(std::move(coords));
@@ -316,7 +326,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
                             //geom::LineString * line = new geom::LineString();
                             if(add_start) {
                                 //line->addPoint(x0,y0);
-                                coords->add(Coordinate(x0, y0));
+                                coords->add(Coordinate(x0, y0, z0));
                                 add_start = false;
                             }
                             //line->addSubLineString(&g, start_index, i-1);
@@ -349,7 +359,7 @@ RectangleIntersection::clip_linestring_parts(const geom::LineString* gi,
                 //geom::LineString * line = new geom::LineString();
                 if(add_start) {
                     //line->addPoint(x0,y0);
-                    coords->add(Coordinate(x0, y0));
+                    coords->add(Coordinate(x0, y0, z0));
                     add_start = false;
                 }
                 //line->addSubLineString(&g, start_index, i-1);
@@ -521,7 +531,6 @@ RectangleIntersection::clip_polygon_to_polygons(const geom::Polygon* g,
 
     parts.reconnectPolygons(rect);
     parts.release(toParts);
-
 }
 
 /**
@@ -682,7 +691,14 @@ std::unique_ptr<geom::Geometry>
 RectangleIntersection::clip(const geom::Geometry& g, const Rectangle& rect)
 {
     RectangleIntersection ri(g, rect);
-    return ri.clip();
+    std::unique_ptr<geom::Geometry> result = ri.clip();
+
+    if (g.hasZ()) {
+        std::unique_ptr<ElevationModel> elevModel = ElevationModel::create(g);
+        elevModel->populateZ(*result);
+    }
+
+    return result;
 }
 
 std::unique_ptr<geom::Geometry>
diff --git a/tests/unit/capi/GEOSClipByRectTest.cpp b/tests/unit/capi/GEOSClipByRectTest.cpp
index 13fde6f732..be47a4c0f7 100644
--- a/tests/unit/capi/GEOSClipByRectTest.cpp
+++ b/tests/unit/capi/GEOSClipByRectTest.cpp
@@ -42,6 +42,27 @@ struct test_capigeosclipbyrect_data : public capitest::utility {
         ensure(equal);
     }
 
+    void
+    checkClipByRectIdentical(
+        const char* wktInput,
+        double xmin, double ymin,
+        double xmax, double ymax,
+        const char* wktExpected)
+    {
+        input_ = fromWKT(wktInput);
+        result_ = GEOSClipByRect(input_, xmin, ymin, xmax, ymax);
+        expected_ = fromWKT(wktExpected);
+
+        bool equal =  GEOSEqualsIdentical(result_, expected_) == 1;
+        if (!equal) {
+            wkt_ = GEOSWKTWriter_write(wktw_, expected_);
+            std::printf("EXP: %s\n", wkt_);
+            GEOSFree(wkt_);
+            wkt_ = GEOSWKTWriter_write(wktw_, result_);
+            std::printf("OBT: %s\n", wkt_);
+        }
+        ensure(equal);
+    }
 };
 
 typedef test_group<test_capigeosclipbyrect_data> group;
@@ -70,7 +91,7 @@ template<>
 template<>
 void object::test<2>()
 {
-    checkClipByRect(
+    checkClipByRectIdentical(
         "POINT(15 15)",
         10, 10, 20, 20,
         "POINT(15 15)"
@@ -106,7 +127,7 @@ template<>
 template<>
 void object::test<5>()
 {
-    checkClipByRect(
+    checkClipByRectIdentical(
         "LINESTRING(15 15, 16 15)",
         10, 10, 20, 20,
         "LINESTRING(15 15, 16 15)"
@@ -130,7 +151,7 @@ template<>
 template<>
 void object::test<7>()
 {
-    checkClipByRect(
+    checkClipByRectIdentical(
         "LINESTRING(10 5, 25 20)",
         10, 10, 20, 20,
         "LINESTRING (15 10, 20 15)"
@@ -191,7 +212,7 @@ template<>
 void object::test<12>()
 {
     const char* wkt = "POLYGON((1 1, 1 30, 30 30, 30 1, 1 1),(10 10, 20 10, 20 20, 10 20, 10 10))";
-    checkClipByRect(
+    checkClipByRectIdentical(
         wkt,
         0, 0, 40, 40,
         wkt);
@@ -202,7 +223,7 @@ template<>
 template<>
 void object::test<13>()
 {
-    checkClipByRect(
+    checkClipByRectIdentical(
         "POLYGON((0 0, 0 30, 30 30, 30 0, 0 0),(10 10, 20 10, 20 20, 10 20, 10 10))",
         5, 5, 15, 15,
         "POLYGON ((5 5, 5 15, 10 15, 10 10, 15 10, 15 5, 5 5))"
@@ -282,6 +303,66 @@ void object::test<17>()
     ensure(result_ == nullptr);
 }
 
+/// Point Z inside
+template<>
+template<>
+void object::test<18>()
+{
+    checkClipByRectIdentical(
+        "POINT Z (15 15 100)",
+        10, 10, 20, 20,
+        "POINT Z (15 15 100)"
+        );
+}
+
+/// Line Z inside
+template<>
+template<>
+void object::test<19>()
+{
+    checkClipByRectIdentical(
+        "LINESTRING Z (15 15 0, 16 15 100)",
+        10, 10, 20, 20,
+        "LINESTRING Z (15 15 0, 16 15 100)"
+        );
+}
+
+/// Line Z splitting rectangle
+template<>
+template<>
+void object::test<20>()
+{
+    checkClipByRectIdentical(
+        "LINESTRING Z (0 15 0, 100 15 100)",
+        10, 10, 20, 20,
+        "LINESTRING Z (10 15 10, 20 15 20)"
+        );
+}
+
+/// Polygon Z overlapping rectangle
+template<>
+template<>
+void object::test<21>()
+{
+    checkClipByRectIdentical(
+        "POLYGON Z ((0 0 100, 0 30 100, 30 30 100, 30 0 100, 0 0 100),(10 10 100, 20 10 100, 20 20 100, 10 20 100, 10 10 100))",
+        5, 5, 15, 15,
+        "POLYGON Z ((5 5 100, 5 15 100, 10 15 100, 10 10 100, 15 10 100, 15 5 100, 5 5 100))"
+        );
+}
+
+/// Polygon Z enclosing rectangle
+template<>
+template<>
+void object::test<22>()
+{
+    checkClipByRectIdentical(
+        "POLYGON Z ((0 0 100, 0 30 100, 30 30 100, 30 0 100, 0 0 100))",
+        5, 5, 15, 15,
+        "POLYGON Z ((5 5 100, 5 15 100, 15 15 100, 15 5 100, 5 5 100))"
+        );
+}
+
 
 } // namespace tut