From 7d9ebd7509080ae91c581678172f6a087c5ba856 Mon Sep 17 00:00:00 2001
From: vadimcn <vadimcn@gmail.com>
Date: Wed, 15 May 2024 07:13:55 -0700
Subject: [PATCH] GEOSDensify: Interpolate Z coordinate (#1094)

---
 include/geos/geom/LineSegment.h     |  3 ++-
 tests/unit/capi/GEOSDensifyTest.cpp | 38 ++++++++++++++++++++++++++++-
 tests/unit/capi/capi_test_utils.h   | 14 +++++++++++
 3 files changed, 53 insertions(+), 2 deletions(-)

diff --git a/include/geos/geom/LineSegment.h b/include/geos/geom/LineSegment.h
index ea1e0a7881..b7f7f712c2 100644
--- a/include/geos/geom/LineSegment.h
+++ b/include/geos/geom/LineSegment.h
@@ -321,7 +321,8 @@ class GEOS_DLL LineSegment {
     {
         ret = Coordinate(
             p0.x + segmentLengthFraction * (p1.x - p0.x),
-            p0.y + segmentLengthFraction * (p1.y - p0.y));
+            p0.y + segmentLengthFraction * (p1.y - p0.y),
+            p0.z + segmentLengthFraction * (p1.z - p0.z));
     };
 
     /** \brief
diff --git a/tests/unit/capi/GEOSDensifyTest.cpp b/tests/unit/capi/GEOSDensifyTest.cpp
index 97ede1605e..721e91dc2c 100644
--- a/tests/unit/capi/GEOSDensifyTest.cpp
+++ b/tests/unit/capi/GEOSDensifyTest.cpp
@@ -29,7 +29,7 @@ struct test_capigeosdensify_data : public capitest::utility
         result_ = GEOSDensify(input_, tolerance);
         ensure("result not NULL", result_ != nullptr);
 
-        ensure_geometry_equals(result_, expected_);
+        ensure_geometry_equals_identical(result_, expected_);
         ensure_equals("result SRID == expected SRID", GEOSGetSRID(result_), srid);
     }
 };
@@ -169,4 +169,40 @@ void object::test<10>()
     result_ = GEOSDensify(input_, 0.1);
     ensure("curved geometries not supported", result_ == nullptr);
 }
+// Densify a LINESTRING Z, check that Z gets interpolated
+template <>
+template <>
+void object::test<11>()
+{
+    testDensify(
+        "LINESTRING Z (0 0 0, 0 6 2)",
+        "LINESTRING Z (0 0 0, 0 3 1, 0 6 2)",
+        3.0
+    );
+}
+
+// Densify a LINEARRING Z
+template <>
+template <>
+void object::test<12>()
+{
+    testDensify(
+        "LINEARRING Z (0 0 0, 0 6 2, 6 6 12, 0 0 0)", 
+        "LINEARRING Z (0 0 0, 0 3 1, 0 6 2, 3 6 7, 6 6 12, 4 4 8, 2 2 4, 0 0 0)",
+        3.0
+    );
+}
+
+// Densify a POLYGON Z
+template <>
+template <>
+void object::test<13>()
+{
+    testDensify(
+        "POLYGON Z ((0 0 0, 10 0 2, 10 10 10, 0 10 2, 0 0 0), (1 1 0, 1 7 0, 7 7 0, 7 1 0, 1 1 0))",
+        "POLYGON Z ((0 0 0, 5 0 1, 10 0 2, 10 5 6, 10 10 10, 5 10 6, 0 10 2, 0 5 1, 0 0 0), (1 1 0, 1 4 0, 1 7 0, 4 7 0, 7 7 0, 7 4 0, 7 1 0, 4 1 0, 1 1 0))",
+        5.0
+    );
+}
+
 } // namespace tut
diff --git a/tests/unit/capi/capi_test_utils.h b/tests/unit/capi/capi_test_utils.h
index ef6f6c07bb..310ed82e34 100644
--- a/tests/unit/capi/capi_test_utils.h
+++ b/tests/unit/capi/capi_test_utils.h
@@ -120,6 +120,20 @@ namespace capitest {
             tut::ensure_equals("GEOSEqualsExact(g1, g2, tolerance)", rslt, 1);
         }
 
+        void
+        ensure_geometry_equals_identical(GEOSGeometry* g1, GEOSGeometry* g2)
+        {
+            char rslt;
+            if (g1 == nullptr || g2 == nullptr) {
+                rslt = (g1 == nullptr && g2 == nullptr) ? 1 : 0;
+            } 
+            else {
+                rslt = GEOSEqualsIdentical(g1, g2);
+            }
+            report_not_equal("ensure_equals_identical", g1, g2, 1e-12, rslt);
+            tut::ensure_equals("GEOSEqualsIdentical(g1, g2)", rslt, 1);
+        }        
+
         void
         ensure_geometry_equals(GEOSGeometry* g1, GEOSGeometry* g2)
         {