Curved geometries: Implement getLength()

libgeos · Jun 20, 2024 · 2647ca0 · 2647ca0
1 parent a6810de
commit 2647ca0
Show file tree

Hide file tree

Showing 11 changed files with 148 additions and 19 deletions.
diff --git a/include/geos/algorithm/CircularArcs.h b/include/geos/algorithm/CircularArcs.h
@@ -23,7 +23,6 @@ namespace algorithm {
 
 class GEOS_DLL CircularArcs {
 public:
-    static constexpr double PI = 3.14159265358979323846;
 
     /// Return the circle center of an arc defined by three points
     static geom::CoordinateXY getCenter(const geom::CoordinateXY& p0, const geom::CoordinateXY& p1,

diff --git a/include/geos/geom/CircularArc.h b/include/geos/geom/CircularArc.h
@@ -47,7 +47,7 @@ class GEOS_DLL CircularArc {
 
     int orientation() const {
         if (!m_orientation_known) {
-            m_orientation = algorithm::Orientation::index(center(), p0, p1);
+            m_orientation = algorithm::Orientation::index(p0, p1, p2);
             m_orientation_known = true;
         }
         return m_orientation;
@@ -71,6 +71,43 @@ class GEOS_DLL CircularArc {
         return m_radius;
     }
 
+    bool isLinear() const {
+        return std::isnan(radius());
+    }
+
+    bool isCircle() const {
+        return p0.equals(p2);
+    }
+
+    double angle() const {
+        if (isCircle()) {
+            return 2*MATH_PI;
+        }
+
+        auto t0 = theta0();
+        auto t2 = theta2();
+
+        if (orientation() == algorithm::Orientation::COUNTERCLOCKWISE) {
+            std::swap(t0, t2);
+        }
+
+        if (t0 < t2) {
+            t0 += 2*MATH_PI;
+        }
+
+        auto diff = t0-t2;
+
+        return diff;
+    }
+
+    double length() const {
+        if (isLinear()) {
+            return p0.distance(p2);
+        }
+
+        return angle()*radius();
+    }
+
     double theta0() const {
         return std::atan2(p0.y - center().y, p0.x - center().x);
     }
@@ -123,11 +160,11 @@ class GEOS_DLL CircularArc {
         t2 -= t0;
         theta -= t0;
 
-        if (t2 < 0){
-            t2 += 2*algorithm::CircularArcs::PI;
+        if (t2 < 0) {
+            t2 += 2*MATH_PI;
         }
         if (theta < 0) {
-            theta += 2*algorithm::CircularArcs::PI;
+            theta += 2*MATH_PI;
         }
 
         return theta >= t2;

diff --git a/include/geos/geom/CircularString.h b/include/geos/geom/CircularString.h
@@ -15,7 +15,6 @@
 #pragma once
 
 #include <geos/geom/SimpleCurve.h>
-#include <geos/util/UnsupportedOperationException.h>
 
 namespace geos {
 namespace geom {
@@ -35,10 +34,7 @@ class GEOS_DLL CircularString : public SimpleCurve {
 
     GeometryTypeId getGeometryTypeId() const override;
 
-    double getLength() const override
-    {
-        throw util::UnsupportedOperationException("Cannot calculate length of CircularString");
-    }
+    double getLength() const override;
 
     bool hasCurvedComponents() const override
     {

diff --git a/src/geom/CircularString.cpp b/src/geom/CircularString.cpp
@@ -12,6 +12,7 @@
  *
  **********************************************************************/
 
+#include <geos/geom/CircularArc.h>
 #include <geos/geom/CircularString.h>
 #include <geos/geom/CoordinateSequence.h>
 #include <geos/geom/GeometryFactory.h>
@@ -49,6 +50,23 @@ CircularString::getGeometryTypeId() const
     return GEOS_CIRCULARSTRING;
 }
 
+double
+CircularString::getLength() const
+{
+    if (isEmpty()) {
+        return 0;
+    }
+
+    const CoordinateSequence& coords = *getCoordinatesRO();
+
+    double tot = 0;
+    for (std::size_t i = 2; i < coords.size(); i += 2) {
+        auto len = CircularArc(coords[i-2], coords[i-1], coords[i]).length();
+        tot += len;
+    }
+    return tot;
+}
+
 CircularString*
 CircularString::reverseImpl() const
 {

diff --git a/tests/unit/capi/GEOSLengthTest.cpp b/tests/unit/capi/GEOSLengthTest.cpp
@@ -1,6 +1,7 @@
 #include <tut/tut.hpp>
 // geos
 #include <geos_c.h>
+#include <geos/constants.h>
 
 #include "capi_test_utils.h"
 
@@ -67,7 +68,8 @@ void object::test<4>()
 
     double length = -1;
     int ret = GEOSLength(input_, &length);
-    ensure_equals("error raised on curved geometry", ret, 0);
+    ensure_equals(ret, 1);
+    ensure_equals(length, geos::MATH_PI);
 }
 
 } // namespace tut

diff --git a/tests/unit/geom/CircularArcTest.cpp b/tests/unit/geom/CircularArcTest.cpp
@@ -0,0 +1,76 @@
+#include <tut/tut.hpp>
+
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/CircularArc.h>
+#include <geos/constants.h>
+
+using geos::geom::CoordinateXY;
+using geos::geom::CircularArc;
+using geos::MATH_PI;
+
+namespace tut {
+
+struct test_circulararc_data {
+
+    const double eps = 1e-8;
+
+    void checkAngle(const CoordinateXY& p0, const CoordinateXY& p1, const CoordinateXY& p2, double expected) {
+        CircularArc arc(p0, p1, p2);
+        ensure_equals(p0.toString() + " / " + p1.toString() + " / " + p2.toString(), arc.angle(), expected, eps);
+
+        CircularArc rev(p2, p1, p0);
+        ensure_equals(p2.toString() + " / " + p1.toString() + " / " + p0.toString(), rev.angle(), expected, eps);
+    }
+
+    void checkLength(const CoordinateXY& p0, const CoordinateXY& p1, const CoordinateXY& p2, double expected) {
+        CircularArc arc(p0, p1, p2);
+        ensure_equals(p0.toString() + " / " + p1.toString() + " / " + p2.toString(), arc.length(), expected, eps);
+
+        CircularArc rev(p2, p1, p0);
+        ensure_equals(p2.toString() + " / " + p1.toString() + " / " + p0.toString(), rev.length(), expected, eps);
+    }
+};
+
+using group = test_group<test_circulararc_data>;
+using object = group::object;
+
+group test_circulararc_group("geos::geom::CircularArc");
+
+// test angle() on unit circle
+template<>
+template<>
+void object::test<1>()
+{
+    auto x = std::sqrt(2.0)/2;
+
+    // full circle
+    checkAngle({-1, 0}, {1, 0}, {-1, 0}, 2*MATH_PI);
+
+    // check half-circles
+    checkAngle({-1, 0}, {0, 1}, {1, 0}, MATH_PI); // top
+    checkAngle({-1, 0}, {0, -1}, {1, 0}, MATH_PI); // bottom
+    checkAngle({0, -1}, {-1, 0}, {0, 1}, MATH_PI); // left
+    checkAngle({0, -1}, {1, 0}, {0, 1}, MATH_PI); // right
+
+    // check quadrants
+    checkAngle({-1, 0}, {-x, x}, {0, 1}, MATH_PI/2); // upper left
+    checkAngle({0, 1}, {x, x}, {1, 0}, MATH_PI/2); // upper right
+    checkAngle({0, -1}, {x, -x}, {1, 0}, MATH_PI/2); // lower right
+    checkAngle({0, -1}, {-x, -x}, {-1, 0}, MATH_PI/2); // lower left
+
+    // 3/4
+    checkAngle({-x, x}, {0, -1}, {x, x}, 1.5*MATH_PI); // mouth up
+    checkAngle({-x, -x}, {0, 1}, {x, -x}, 1.5*MATH_PI); // mouth down
+    checkAngle({-x, x}, {1, 0}, {-x, -x}, 1.5*MATH_PI); // mouth left
+    checkAngle({x, -x}, {-1, 0}, {x, x}, 1.5*MATH_PI); // mouth right
+}
+
+// test length()
+template<>
+template<>
+void object::test<2>()
+{
+    checkLength({1.6, 0.4}, {1.6, 0.5}, {1.7, 1}, 0.6122445326877711);
+}
+
+}
diff --git a/tests/unit/geom/CircularStringTest.cpp b/tests/unit/geom/CircularStringTest.cpp
@@ -58,7 +58,7 @@ void object::test<1>()
     ensure(cs->getCoordinate() == nullptr);
 
     ensure_equals(cs->getArea(), 0);
-    ensure_THROW(cs_->getLength(), geos::util::UnsupportedOperationException);
+    ensure_equals(cs->getLength(), 0);
 }
 
 // Basic Geometry API
@@ -74,7 +74,7 @@ void object::test<2>()
     // Geometry size functions
     ensure("isEmpty", !cs_->isEmpty());
     ensure_equals("getArea", cs_->getArea(), 0);
-    ensure_THROW(cs_->getLength(), geos::util::UnsupportedOperationException);
+    ensure_equals("getLength", cs_->getLength(), 2*geos::MATH_PI);
     ensure_equals("getNumGeometries", cs_->getNumGeometries(), 1u);
     ensure_equals("getNumPoints", cs_->getNumPoints(), 5u);
     geos::geom::Envelope expected(0, 4, -1, 1);

diff --git a/tests/unit/geom/CompoundCurveTest.cpp b/tests/unit/geom/CompoundCurveTest.cpp
@@ -65,7 +65,7 @@ void object::test<1>()
     ensure("getCoordinate", cc->getCoordinate() == nullptr);
 
     ensure_equals("getArea", cc->getArea(), 0);
-    ensure_THROW(cc_->getLength(), geos::util::UnsupportedOperationException);
+    ensure_equals("getLength", cc->getLength(), 0);
 }
 
 // Basic Geometry API
@@ -81,7 +81,7 @@ void object::test<2>()
     // Geometry size functions
     ensure("isEmpty", !cc_->isEmpty());
     ensure_equals("getArea", cc_->getArea(), 0);
-    ensure_THROW(cc_->getLength(), geos::util::UnsupportedOperationException);
+    ensure_equals("getLength", cc_->getLength(), geos::MATH_PI + 2);
     ensure_equals("getNumGeometries", cc_->getNumGeometries(), 1u);
     ensure_equals("getNumCurves", cc_->getNumCurves(), 2u);
     ensure_equals("getNumPoints", cc_->getNumPoints(), 5u); // FIXME should this be 5 or 4?

diff --git a/tests/unit/geom/CurvePolygonTest.cpp b/tests/unit/geom/CurvePolygonTest.cpp
@@ -91,7 +91,7 @@ void object::test<2>()
     // Geometry size functions
     ensure("isEmpty", !cp_->isEmpty());
     ensure_THROW(cp_->getArea(), geos::util::UnsupportedOperationException);
-    ensure_THROW(cp_->getLength(), geos::util::UnsupportedOperationException);
+    ensure_equals("getLength", cp_->getLength(), 19.236489581872586, 1e-8);
     ensure_equals("getNumGeometries", cp_->getNumGeometries(), 1u);
     ensure_equals("getNumPoints", cp_->getNumPoints(), 14u);
     ensure_equals("getNumInteriorRing", cp_->getNumInteriorRing(), 1u);

diff --git a/tests/unit/geom/MultiCurveTest.cpp b/tests/unit/geom/MultiCurveTest.cpp
@@ -81,7 +81,7 @@ void object::test<1>()
     ensure("getCoordinate", mc->getCoordinate() == nullptr);
 
     ensure_equals("getArea", mc->getArea(), 0);
-    ensure_THROW(mc_->getLength(), geos::util::UnsupportedOperationException);
+    ensure_equals("getLength", mc->getLength(), 0);
 }
 
 // Basic Geometry API
@@ -97,7 +97,7 @@ void object::test<2>()
     // Geometry size functions
     ensure("isEmpty", !mc_->isEmpty());
     ensure_equals("getArea", mc_->getArea(), 0);
-    ensure_THROW(mc_->getLength(), geos::util::UnsupportedOperationException);
+    ensure_equals("getLength", mc_->getLength(), 22.064916706618778, 1e-8);
     ensure_equals("getNumGeometries", mc_->getNumGeometries(), 3u);
     ensure_equals("getNumPoints", mc_->getNumPoints(), 16u);
     ensure(!mc_->getEnvelopeInternal()->isNull());

diff --git a/tests/unit/geom/MultiSurfaceTest.cpp b/tests/unit/geom/MultiSurfaceTest.cpp
@@ -7,6 +7,7 @@
 #include <geos/geom/MultiSurface.h>
 #include <geos/geom/IntersectionMatrix.h>
 #include <geos/io/WKTReader.h>
+#include <geos/util/UnsupportedOperationException.h>
 
 using geos::geom::CoordinateXY;
 
@@ -85,7 +86,7 @@ void object::test<2>()
     // Geometry size functions
     ensure("isEmpty", !ms_->isEmpty());
     ensure_THROW(ms_->getArea(), geos::util::UnsupportedOperationException);
-    ensure_THROW(ms_->getLength(), geos::util::UnsupportedOperationException);
+    ensure_equals("getLength", ms_->getLength(), 10.283185307179586);
     ensure_equals("getNumGeometries", ms_->getNumGeometries(), 2u);
     ensure_equals("getNumPoints", ms_->getNumPoints(), 10u);
     ensure(!ms_->getEnvelopeInternal()->isNull());