From f7f47c509d9cf4279bb532099cd693d4502a8a2b Mon Sep 17 00:00:00 2001 From: Daniel Baston Date: Thu, 20 Jun 2024 17:31:24 -0400 Subject: [PATCH] Curved geometries: Implement getArea --- include/geos/algorithm/Area.h | 7 ++++ include/geos/geom/CircularArc.h | 10 +++++ include/geos/geom/CircularString.h | 4 ++ include/geos/geom/LineString.h | 4 ++ include/geos/geom/SimpleCurve.h | 2 + src/algorithm/Area.cpp | 55 ++++++++++++++++++++++++ src/geom/CurvePolygon.cpp | 7 +++- tests/unit/algorithm/AreaTest.cpp | 62 +++++++++++++++++++++++----- tests/unit/geom/CircularArcTest.cpp | 15 +++++++ tests/unit/geom/CurvePolygonTest.cpp | 4 +- tests/unit/geom/MultiSurfaceTest.cpp | 2 +- 11 files changed, 157 insertions(+), 15 deletions(-) diff --git a/include/geos/algorithm/Area.h b/include/geos/algorithm/Area.h index f53938c1af..c3c799f877 100644 --- a/include/geos/algorithm/Area.h +++ b/include/geos/algorithm/Area.h @@ -23,12 +23,19 @@ #include namespace geos { + +namespace geom { +class Curve; +} + namespace algorithm { // geos::algorithm class GEOS_DLL Area { public: + static double ofClosedCurve(const geom::Curve& ring); + /** * Computes the area for a ring. * diff --git a/include/geos/geom/CircularArc.h b/include/geos/geom/CircularArc.h index 69ca667460..1a87583a36 100644 --- a/include/geos/geom/CircularArc.h +++ b/include/geos/geom/CircularArc.h @@ -108,6 +108,16 @@ class GEOS_DLL CircularArc { return angle()*radius(); } + double getArea() const { + if (isLinear()) { + return 0; + } + + auto R = radius(); + auto theta = angle(); + return R*R/2*(theta - std::sin(theta)); + } + double theta0() const { return std::atan2(p0.y - center().y, p0.x - center().x); } diff --git a/include/geos/geom/CircularString.h b/include/geos/geom/CircularString.h index fd43bc2b9f..85dc3014e3 100644 --- a/include/geos/geom/CircularString.h +++ b/include/geos/geom/CircularString.h @@ -41,6 +41,10 @@ class GEOS_DLL CircularString : public SimpleCurve { return true; } + bool isCurved() const override { + return true; + } + std::unique_ptr reverse() const { return std::unique_ptr(reverseImpl()); diff --git a/include/geos/geom/LineString.h b/include/geos/geom/LineString.h index faf021f205..478970874a 100644 --- a/include/geos/geom/LineString.h +++ b/include/geos/geom/LineString.h @@ -92,6 +92,10 @@ class GEOS_DLL LineString: public SimpleCurve { double getLength() const override; + bool isCurved() const override { + return false; + } + /** * Creates a LineString whose coordinates are in the reverse * order of this object's diff --git a/include/geos/geom/SimpleCurve.h b/include/geos/geom/SimpleCurve.h index a41ea0ee55..1bf674f694 100644 --- a/include/geos/geom/SimpleCurve.h +++ b/include/geos/geom/SimpleCurve.h @@ -94,6 +94,8 @@ class GEOS_DLL SimpleCurve : public Curve { virtual bool isCoordinate(CoordinateXY& pt) const; + virtual bool isCurved() const = 0; + bool isEmpty() const override; /** \brief diff --git a/src/algorithm/Area.cpp b/src/algorithm/Area.cpp index e51954dbcf..3d2375c957 100644 --- a/src/algorithm/Area.cpp +++ b/src/algorithm/Area.cpp @@ -20,6 +20,10 @@ #include #include +#include +#include +#include +#include using geos::geom::CoordinateXY; @@ -93,7 +97,58 @@ Area::ofRingSigned(const geom::CoordinateSequence* ring) return sum / 2.0; } +double +Area::ofClosedCurve(const geom::Curve& ring) { + if (!ring.isClosed()) { + throw util::IllegalArgumentException("Argument is not closed"); + } + + double sum = 0; + + for (std::size_t i = 0; i < ring.getNumCurves(); i++) { + const geom::SimpleCurve& section = *ring.getCurveN(i); + + if (section.isEmpty()) { + continue; + } + + const geom::CoordinateSequence& coords = *section.getCoordinatesRO(); + + if (section.isCurved()) { + for (std::size_t j = 2; j < coords.size(); j += 2) { + const CoordinateXY& p0 = coords.getAt(j-2); + const CoordinateXY& p1 = coords.getAt(j-1); + const CoordinateXY& p2 = coords.getAt(j); + double triangleArea = 0.5*(p0.x*p2.y - p2.x*p0.y); + sum += triangleArea; + + geom::CircularArc arc(p0, p1, p2); + if (arc.isLinear()) { + continue; + } + + double circularSegmentArea = arc.getArea(); + + if (algorithm::Orientation::index(p0, p2, p1) == algorithm::Orientation::CLOCKWISE) { + sum += circularSegmentArea; + } else { + sum -= circularSegmentArea; + } + } + } else { + for (std::size_t j = 1; j < coords.size(); j++) { + const CoordinateXY& p0 = coords.getAt(j-1); + const CoordinateXY& p1 = coords.getAt(j); + + double triangleArea = 0.5*(p0.x*p1.y - p1.x*p0.y); + sum += triangleArea; + } + } + } + + return std::abs(sum); +} } // namespace geos.algorithm } //namespace geos diff --git a/src/geom/CurvePolygon.cpp b/src/geom/CurvePolygon.cpp index 74a5e943b7..c5213602c7 100644 --- a/src/geom/CurvePolygon.cpp +++ b/src/geom/CurvePolygon.cpp @@ -12,6 +12,7 @@ * **********************************************************************/ +#include #include #include #include @@ -53,7 +54,11 @@ namespace geom { } double CurvePolygon::getArea() const { - throw util::UnsupportedOperationException(); + double sum = algorithm::Area::ofClosedCurve(*shell); + for (const auto& hole : holes) { + sum -= algorithm::Area::ofClosedCurve(*hole); + } + return sum; } bool CurvePolygon::hasCurvedComponents() const { diff --git a/tests/unit/algorithm/AreaTest.cpp b/tests/unit/algorithm/AreaTest.cpp index 1a81ff08b9..2033f85b1a 100644 --- a/tests/unit/algorithm/AreaTest.cpp +++ b/tests/unit/algorithm/AreaTest.cpp @@ -19,6 +19,7 @@ using namespace geos; using namespace geos::geom; +using geos::algorithm::Area; namespace tut { // @@ -33,7 +34,7 @@ struct test_area_data { geos::io::WKTReader reader_; test_area_data(): geom_(nullptr), - pm_(1), + pm_(), factory_(GeometryFactory::create(&pm_, 0)), reader_(factory_.get()) { assert(nullptr == geom_); @@ -48,19 +49,23 @@ struct test_area_data { void checkAreaOfRing(std::string wkt, double expectedArea) { - std::unique_ptr lineGeom(reader_.read(wkt)); - std::unique_ptr line(dynamic_cast(lineGeom.release())); - ensure(nullptr != line.get()); - const CoordinateSequence* ringSeq = line->getCoordinatesRO(); + auto ringGeom = reader_.read(wkt); - std::vector ringCoords; - ringSeq->toVector(ringCoords); + if (const LineString* line = dynamic_cast(ringGeom.get())) { + const CoordinateSequence* ringSeq = line->getCoordinatesRO(); - double actual1 = algorithm::Area::ofRing(ringCoords); - double actual2 = algorithm::Area::ofRing(ringSeq); + std::vector ringCoords; + ringSeq->toVector(ringCoords); - ensure_equals(actual1, expectedArea); - ensure_equals(actual2, expectedArea); + double actual1 = algorithm::Area::ofRing(ringCoords); + double actual2 = algorithm::Area::ofRing(ringSeq); + + ensure_equals(actual1, expectedArea); + ensure_equals(actual2, expectedArea); + } + + double actual3 = algorithm::Area::ofClosedCurve(*ringGeom); + ensure_equals("Area::ofClosedCurve", actual3, expectedArea, 1e-6); } void @@ -115,6 +120,41 @@ void object::test<3> checkAreaOfRingSigned("LINESTRING (100 200, 100 100, 200 100, 200 200, 100 200)", -10000.0); } +template<> +template<> +void object::test<4> +() +{ + checkAreaOfRing("CIRCULARSTRING (0 0, 2 2, 4 0, 2 -2, 0 0)", 4*MATH_PI); +} + +template<> +template<> +void object::test<5> +() +{ + checkAreaOfRing("COMPOUNDCURVE (CIRCULARSTRING (0 0, 2 2, 4 0), (4 0, 0 0))", 2*MATH_PI); +} + +template<> +template<> +void object::test<6> +() +{ + // expected ares from PostGIS after ST_CurveToLine(geom, 1e-13, 1) + checkAreaOfRing("CIRCULARSTRING (0 0, 2 2, 4 0, 2 1, 0 0)", 3.48759); +} + +template<> +template<> +void object::test<7> +() +{ + // expected ares from PostGIS after ST_CurveToLine(geom, 1e-13, 1) + checkAreaOfRing("COMPOUNDCURVE (CIRCULARSTRING (0 0, 2 0, 2 1, 2 3, 4 3, 4 5, 1 4, 0.5 0.8, 0 0))", 11.243342); + checkAreaOfRing("COMPOUNDCURVE (CIRCULARSTRING (0 0, 2 0, 2 1, 2 3, 4 3), (4 3, 4 5, 1 4, 0 0))", 9.321903); +} + } // namespace tut diff --git a/tests/unit/geom/CircularArcTest.cpp b/tests/unit/geom/CircularArcTest.cpp index 00f7bbcbd3..a5d3c57726 100644 --- a/tests/unit/geom/CircularArcTest.cpp +++ b/tests/unit/geom/CircularArcTest.cpp @@ -73,4 +73,19 @@ void object::test<2>() checkLength({1.6, 0.4}, {1.6, 0.5}, {1.7, 1}, 0.6122445326877711); } + +// test getArea() +template<> +template<> +void object::test<3>() +{ + ensure_equals("half circle, R=2", CircularArc({-2, 0}, {0, 2}, {2, 0}).getArea(), MATH_PI*2); + + ensure_equals("full circle, R=3", CircularArc({-3, 0}, {3, 0}, {-3, 0}).getArea(), MATH_PI*3*3); + + ensure_equals("3/4, mouth up, R=2", CircularArc({-std::sqrt(2), std::sqrt(2)}, {0, -2}, {std::sqrt(2), std::sqrt(2)}).getArea(), MATH_PI*4 - 2*(MATH_PI/2-1), 1e-8); + + ensure_equals("1/4, pointing up, R=2", CircularArc({-std::sqrt(2), std::sqrt(2)}, {0, 2}, {std::sqrt(2), std::sqrt(2)}).getArea(), 2*(MATH_PI/2-1), 1e-8); +} + } diff --git a/tests/unit/geom/CurvePolygonTest.cpp b/tests/unit/geom/CurvePolygonTest.cpp index 9cafc6516b..0ee7d4baa6 100644 --- a/tests/unit/geom/CurvePolygonTest.cpp +++ b/tests/unit/geom/CurvePolygonTest.cpp @@ -74,7 +74,7 @@ void object::test<1>() ensure("getCoordinates", cp->getCoordinates()->isEmpty()); ensure("getCoordinate", cp->getCoordinate() == nullptr); - ensure_THROW(cp->getArea(), geos::util::UnsupportedOperationException); + ensure_equals("getArea", cp->getArea(), 0.0); ensure_equals("getLength", cp->getLength(), 0.0); } @@ -90,7 +90,7 @@ void object::test<2>() // Geometry size functions ensure("isEmpty", !cp_->isEmpty()); - ensure_THROW(cp_->getArea(), geos::util::UnsupportedOperationException); + ensure_equals("getArea", cp_->getArea(), 9.0526564962674, 1e-8); // expected value from PostGIS with ST_CurveToLine(geom, 1e-13, 1) ensure_equals("getLength", cp_->getLength(), 19.236489581872586, 1e-8); ensure_equals("getNumGeometries", cp_->getNumGeometries(), 1u); ensure_equals("getNumPoints", cp_->getNumPoints(), 14u); diff --git a/tests/unit/geom/MultiSurfaceTest.cpp b/tests/unit/geom/MultiSurfaceTest.cpp index 751010c23c..1492f22792 100644 --- a/tests/unit/geom/MultiSurfaceTest.cpp +++ b/tests/unit/geom/MultiSurfaceTest.cpp @@ -85,7 +85,7 @@ void object::test<2>() // Geometry size functions ensure("isEmpty", !ms_->isEmpty()); - ensure_THROW(ms_->getArea(), geos::util::UnsupportedOperationException); + ensure_equals("getArea", ms_->getArea(), 4.141592653589132, 1e-6); // expected value from PostGIS with ST_CurveToLine(geom, 1e-13, 1) ensure_equals("getLength", ms_->getLength(), 10.283185307179586); ensure_equals("getNumGeometries", ms_->getNumGeometries(), 2u); ensure_equals("getNumPoints", ms_->getNumPoints(), 10u);