From d9bee63471dea5596c5aaa89eaa0a8512f937773 Mon Sep 17 00:00:00 2001 From: Daniel Baston Date: Fri, 12 Apr 2024 12:49:25 -0400 Subject: [PATCH] Implement getEnvelopeInternal for curved types --- include/geos/algorithm/CircularArcs.h | 37 ++++ include/geos/geom/CircularString.h | 13 ++ include/geos/geom/LineString.h | 5 + include/geos/geom/SimpleCurve.h | 8 +- src/algorithm/CircularArcs.cpp | 124 ++++++++++++ src/geom/CircularString.cpp | 9 + src/geom/LineString.cpp | 2 +- src/geom/SimpleCurve.cpp | 19 +- tests/unit/algorithm/CircularArcsTest.cpp | 219 ++++++++++++++++++++++ tests/unit/geom/CircularStringTest.cpp | 4 +- tests/unit/geom/CompoundCurveTest.cpp | 3 +- tests/unit/geom/CurvePolygonTest.cpp | 10 +- 12 files changed, 439 insertions(+), 14 deletions(-) create mode 100644 include/geos/algorithm/CircularArcs.h create mode 100644 src/algorithm/CircularArcs.cpp create mode 100644 tests/unit/algorithm/CircularArcsTest.cpp diff --git a/include/geos/algorithm/CircularArcs.h b/include/geos/algorithm/CircularArcs.h new file mode 100644 index 0000000000..22293e5482 --- /dev/null +++ b/include/geos/algorithm/CircularArcs.h @@ -0,0 +1,37 @@ +/********************************************************************** + * + * GEOS - Geometry Engine Open Source + * http://geos.osgeo.org + * + * Copyright (C) 2024 ISciences, LLC + * + * This is free software; you can redistribute and/or modify it under + * the terms of the GNU Lesser General Public Licence as published + * by the Free Software Foundation. + * See the COPYING file for more information. + * + **********************************************************************/ + +#pragma once + +#include +#include + +#include "geos/geom/Envelope.h" + +namespace geos { +namespace algorithm { + +class GEOS_DLL CircularArcs { +public: + /// Return the circle center of an arc defined by three points + static geom::CoordinateXY getCenter(const geom::CoordinateXY& p0, const geom::CoordinateXY& p1, + const geom::CoordinateXY& p2); + + /// Expand an envelope to include an arc defined by three points + static void expandEnvelope(geom::Envelope& e, const geom::CoordinateXY& p0, const geom::CoordinateXY& p1, + const geom::CoordinateXY& p2); +}; + +} +} diff --git a/include/geos/geom/CircularString.h b/include/geos/geom/CircularString.h index e600092250..7d34cd9b2e 100644 --- a/include/geos/geom/CircularString.h +++ b/include/geos/geom/CircularString.h @@ -56,6 +56,19 @@ class GEOS_DLL CircularString : public SimpleCurve { return true; } +protected: + + /// \brief + /// Constructs a CircularString taking ownership the + /// given CoordinateSequence. + CircularString(std::unique_ptr && pts, + const GeometryFactory& newFactory); + + void geometryChangedAction() override + { + envelope = computeEnvelopeInternal(false); + } + }; diff --git a/include/geos/geom/LineString.h b/include/geos/geom/LineString.h index ebcaffc3da..faf021f205 100644 --- a/include/geos/geom/LineString.h +++ b/include/geos/geom/LineString.h @@ -120,6 +120,11 @@ class GEOS_DLL LineString: public SimpleCurve { return SORTINDEX_LINESTRING; }; + void geometryChangedAction() override + { + envelope = computeEnvelopeInternal(true); + } + private: void validateConstruction(); diff --git a/include/geos/geom/SimpleCurve.h b/include/geos/geom/SimpleCurve.h index 23c9467429..187af0c91e 100644 --- a/include/geos/geom/SimpleCurve.h +++ b/include/geos/geom/SimpleCurve.h @@ -126,14 +126,10 @@ class GEOS_DLL SimpleCurve : public Curve { SimpleCurve(const SimpleCurve& other); SimpleCurve(std::unique_ptr&& newCoords, + bool isLinear, const GeometryFactory& factory); - Envelope computeEnvelopeInternal() const; - - void geometryChangedAction() override - { - envelope = computeEnvelopeInternal(); - } + Envelope computeEnvelopeInternal(bool isLinear) const; // TODO: hold value or shared_ptr instead of unique_ptr std::unique_ptr points; diff --git a/src/algorithm/CircularArcs.cpp b/src/algorithm/CircularArcs.cpp new file mode 100644 index 0000000000..abfd73f42c --- /dev/null +++ b/src/algorithm/CircularArcs.cpp @@ -0,0 +1,124 @@ +/********************************************************************** + * + * GEOS - Geometry Engine Open Source + * http://geos.osgeo.org + * + * Copyright (C) 2024 ISciences, LLC + * + * This is free software; you can redistribute and/or modify it under + * the terms of the GNU Lesser General Public Licence as published + * by the Free Software Foundation. + * See the COPYING file for more information. + * + **********************************************************************/ + +#include + +#include "geos/algorithm/Angle.h" +#include "geos/geom/Envelope.h" +#include "geos/geom/Quadrant.h" + +using geos::geom::CoordinateXY; + +namespace geos { +namespace algorithm { + +CoordinateXY +CircularArcs::getCenter(const CoordinateXY& p0, const CoordinateXY& p1, const CoordinateXY& p2) +{ + if (p0.equals2D(p2)) { + // Closed circle + return { 0.5*(p0.x + p1.x), 0.5*(p0.y + p1.y) }; + } + + // Circumcenter formulas from Graphics Gems III + CoordinateXY a{p1.x - p2.x, p1.y - p2.y}; + CoordinateXY b{p2.x - p0.x, p2.y - p0.y}; + CoordinateXY c{p0.x - p1.x, p0.y - p1.y}; + + double d1 = -(b.x*c.x + b.y*c.y); + double d2 = -(c.x*a.x + c.y*a.y); + double d3 = -(a.x*b.x + a.y*b.y); + + double e1 = d2*d3; + double e2 = d3*d1; + double e3 = d1*d2; + double e = e1 + e2 + e3; + + CoordinateXY G3{p0.x + p1.x + p2.x, p0.y + p1.y + p2.y}; + CoordinateXY H {(e1*p0.x + e2*p1.x + e3*p2.x) / e, (e1*p0.y + e2*p1.y + e3*p2.y) / e}; + + CoordinateXY center = {0.5*(G3.x - H.x), 0.5*(G3.y - H.y)}; + + return center; +} + +void +CircularArcs::expandEnvelope(geom::Envelope& e, const geom::CoordinateXY& p0, const geom::CoordinateXY& p1, + const geom::CoordinateXY& p2) +{ + e.expandToInclude(p0); + e.expandToInclude(p1); + e.expandToInclude(p2); + + CoordinateXY center = getCenter(p0, p1, p2); + + // zero-length arc + if (center.equals2D(p0) || center.equals2D(p1)) { + return; + } + + // collinear + if (std::isnan(center.x)) { + return; + } + + auto orientation = Orientation::index(center, p0, p1); + + //* 1 | 0 + //* --+-- + //* 2 | 3 + + using geom::Quadrant; + + auto q0 = geom::Quadrant::quadrant(center, p0); + auto q2 = geom::Quadrant::quadrant(center, p2); + double R = center.distance(p1); + + if (q0 == q2) { + // Start and end quadrants are the same. Either the arc crosses all of + // the axes, or none of the axes. + if (Orientation::index(center, p1, p2) != orientation) { + e.expandToInclude({center.x, center.y + R}); + e.expandToInclude({center.x - R, center.y}); + e.expandToInclude({center.x, center.y - R}); + e.expandToInclude({center.x + R, center.y}); + } + + return; + } + + if (orientation == Orientation::CLOCKWISE) { + std::swap(q0, q2); + } + + for (auto q = q0 + 1; (q % 4) != ((q2+1) % 4); q++) { + switch (q % 4) { + case Quadrant::NW: + e.expandToInclude({center.x, center.y + R}); + break; + case Quadrant::SW: + e.expandToInclude({center.x - R, center.y}); + break; + case Quadrant::SE: + e.expandToInclude({center.x, center.y - R}); + break; + case Quadrant::NE: + e.expandToInclude({center.x + R, center.y}); + break; + } + } +} + +} +} diff --git a/src/geom/CircularString.cpp b/src/geom/CircularString.cpp index a4376fcb36..1bdaaf16d8 100644 --- a/src/geom/CircularString.cpp +++ b/src/geom/CircularString.cpp @@ -20,6 +20,15 @@ namespace geos { namespace geom { +/*public*/ +CircularString::CircularString(std::unique_ptr && newCoords, + const GeometryFactory& factory) + : + SimpleCurve(std::move(newCoords), false, factory) +{ + //validateConstruction(); +} + CircularString::~CircularString() = default; std::string diff --git a/src/geom/LineString.cpp b/src/geom/LineString.cpp index f3a97822b0..90e37479f6 100644 --- a/src/geom/LineString.cpp +++ b/src/geom/LineString.cpp @@ -58,7 +58,7 @@ LineString::LineString(const LineString& ls) LineString::LineString(CoordinateSequence::Ptr && newCoords, const GeometryFactory& factory) : - SimpleCurve(std::move(newCoords), factory) + SimpleCurve(std::move(newCoords), true, factory) { validateConstruction(); } diff --git a/src/geom/SimpleCurve.cpp b/src/geom/SimpleCurve.cpp index 68c8414357..a0e0397608 100644 --- a/src/geom/SimpleCurve.cpp +++ b/src/geom/SimpleCurve.cpp @@ -17,6 +17,7 @@ #include +#include #include #include #include @@ -38,21 +39,33 @@ SimpleCurve::SimpleCurve(const SimpleCurve& other) } SimpleCurve::SimpleCurve(std::unique_ptr&& newCoords, + bool isLinear, const GeometryFactory& factory) : Curve(factory), points(newCoords ? std::move(newCoords) : std::make_unique()), - envelope(computeEnvelopeInternal()) + envelope(computeEnvelopeInternal(isLinear)) { } Envelope -SimpleCurve::computeEnvelopeInternal() const +SimpleCurve::computeEnvelopeInternal(bool isLinear) const { if(isEmpty()) { return Envelope(); } - return points->getEnvelope(); + if(isLinear) { + return points->getEnvelope(); + } else { + Envelope e; + for (std::size_t i = 2; i < points->size(); i++) { + algorithm::CircularArcs::expandEnvelope(e, + points->getAt(i-2), + points->getAt(i-1), + points->getAt(i)); + } + return e; + } } std::unique_ptr diff --git a/tests/unit/algorithm/CircularArcsTest.cpp b/tests/unit/algorithm/CircularArcsTest.cpp new file mode 100644 index 0000000000..475261d4fb --- /dev/null +++ b/tests/unit/algorithm/CircularArcsTest.cpp @@ -0,0 +1,219 @@ +#include + +#include +#include + +using geos::geom::CoordinateXY; +using geos::algorithm::CircularArcs; +using geos::geom::Envelope; + +namespace tut { + +struct test_circulararcs_data { + const double eps = 1e-8; + + void checkEnvelope(const CoordinateXY& p0, const CoordinateXY& p1, const CoordinateXY& p2, + double xmin, double ymin, double xmax, double ymax) + { + { + Envelope e; + CircularArcs::expandEnvelope(e, p0, p1, p2); + + ensure_equals("p0-p1-p2 xmin", e.getMinX(), xmin, eps); + ensure_equals("p0-p1-p2 xmax", e.getMaxX(), xmax, eps); + ensure_equals("p0-p1-p2 ymin", e.getMinY(), ymin, eps); + ensure_equals("p0-p1-p2 ymax", e.getMaxY(), ymax, eps); + } + + { + Envelope e; + CircularArcs::expandEnvelope(e, p2, p1, p0); + + ensure_equals("p2-p1-p0 xmin", e.getMinX(), xmin, eps); + ensure_equals("p2-p1-p0 xmax", e.getMaxX(), xmax, eps); + ensure_equals("p2-p1-p0 ymin", e.getMinY(), ymin, eps); + ensure_equals("p2-p1-p0 ymax", e.getMaxY(), ymax, eps); + } + } +}; + +using group = test_group; +using object = group::object; + +group test_circulararcs_group("geos::algorithm::CircularArcs"); + +template<> +template<> +void object::test<1>() +{ + CoordinateXY p0{0, 10}; + CoordinateXY p1{100, 110}; + CoordinateXY p2{200, 10}; + + auto center = CircularArcs::getCenter(p0, p1, p2); + + ensure_equals(center, CoordinateXY{100, 10}); +} + +template<> +template<> +void object::test<2>() +{ + CoordinateXY p0{0, 0}; + CoordinateXY p1{1, 1}; + CoordinateXY p2{0, 2}; + + auto center = CircularArcs::getCenter(p0, p1, p2); + + ensure_equals(center, CoordinateXY{0, 1}); +} + +template<> +template<> +void object::test<3>() +{ + CoordinateXY p0{54.22, 31.8}; + CoordinateXY p1{16.07, 11.9}; + CoordinateXY p2{12.22, 3.99}; + + auto center = CircularArcs::getCenter(p0, p1, p2); + + ensure(center.distance(CoordinateXY{52.0123, -10.486}) < 1e-4); +} + +// complete circle +template<> +template<> +void object::test<4>() +{ + CoordinateXY p0{3, 4}; + CoordinateXY p1{7, 8}; + CoordinateXY p2{3, 4}; + + auto center = CircularArcs::getCenter(p0, p1, p2); + + ensure_equals(center, CoordinateXY{5, 6}); +} + +// collinear +template<> +template<> +void object::test<5>() +{ + CoordinateXY p0{1, 2}; + CoordinateXY p1{2, 3}; + CoordinateXY p2{3, 4}; + + auto center = CircularArcs::getCenter(p0, p1, p2); + + ensure(std::isnan(center.x)); + ensure(std::isnan(center.y)); +} + +// CCW quadrant 2 to quadrant 1 +template<> +template<> +void object::test<6>() +{ + CoordinateXY p0{-std::sqrt(2), -std::sqrt(2)}; + CoordinateXY p1{2, 0}; + CoordinateXY p2{-std::sqrt(2), std::sqrt(2)}; + + checkEnvelope(p0, p1, p2, + -std::sqrt(2), -2, 2, 2); +} + +// quadrant 0 to quadrant 0, crossing all axes +template<> +template<> +void object::test<7>() +{ + CoordinateXY p0{std::sqrt(2), std::sqrt(2)}; + CoordinateXY p1{2, 0}; + CoordinateXY p2{std::sqrt(3), 1}; + + checkEnvelope(p0, p1, p2, -2, -2, 2, 2); +} + + +// quadrant 0 to quadrant 0, crossing no axes +template<> +template<> +void object::test<8>() +{ + CoordinateXY p0{1, std::sqrt(3)}; + CoordinateXY p1{std::sqrt(2), std::sqrt(2)}; + CoordinateXY p2{std::sqrt(3), 1}; + + checkEnvelope(p0, p1, p2, + 1, 1, std::sqrt(3), std::sqrt(3)); +} + +// half circle with start points on -/+ x axis +template<> +template<> +void object::test<9>() +{ + CoordinateXY p0{-1, 0}; + CoordinateXY p1{0, 1}; + CoordinateXY p2{1, 0}; + + checkEnvelope(p0, p1, p2, + -1, 0, 1, 1); +} + +// CCW quadrant 0 to quadrant 3 +template<> +template<> +void object::test<10>() +{ + CoordinateXY p0{std::sqrt(2), std::sqrt(2)}; + CoordinateXY p1{-2, 0}; + CoordinateXY p2{std::sqrt(2), -std::sqrt(2)}; + + checkEnvelope(p0, p1, p2, + -2, -2, std::sqrt(2), 2); +} + +// collinear +template<> +template<> +void object::test<11>() +{ + CoordinateXY p0{-1, -1}; + CoordinateXY p1{1, 1}; + CoordinateXY p2{2, 2}; + + checkEnvelope(p0, p1, p2, + -1, -1, 2, 2); +} + +// collinear +template<> +template<> +void object::test<12>() +{ + + CoordinateXY p0{1, 2}; + CoordinateXY p1{2, 3}; + CoordinateXY p2{3, 4}; + + checkEnvelope(p0, p1, p2, + 1, 2, 3, 4); +} + +// repeated +template<> +template<> +void object::test<13>() +{ + CoordinateXY p0{3, 4}; + CoordinateXY p1{3, 4}; + CoordinateXY p2{3, 4}; + + checkEnvelope(p0, p1, p2, + 3, 4, 3, 4); +} + +} + diff --git a/tests/unit/geom/CircularStringTest.cpp b/tests/unit/geom/CircularStringTest.cpp index 6f5cf03ae3..8e1eaf4f3e 100644 --- a/tests/unit/geom/CircularStringTest.cpp +++ b/tests/unit/geom/CircularStringTest.cpp @@ -77,8 +77,8 @@ void object::test<2>() ensure_THROW(cs_->getLength(), geos::util::UnsupportedOperationException); ensure_equals("getNumGeometries", cs_->getNumGeometries(), 1u); ensure_equals("getNumPoints", cs_->getNumPoints(), 5u); - ensure(!cs_->getEnvelopeInternal()->isNull()); - // FIXME calculate envelope correctly + geos::geom::Envelope expected(0, 4, -1, 1); + ensure("getEnvelopeInternal", cs_->getEnvelopeInternal()->equals(&expected)); // Geometry dimension functions ensure_equals("getDimension", cs_->getDimension(), geos::geom::Dimension::L); diff --git a/tests/unit/geom/CompoundCurveTest.cpp b/tests/unit/geom/CompoundCurveTest.cpp index c394bf945c..ae232f5caa 100644 --- a/tests/unit/geom/CompoundCurveTest.cpp +++ b/tests/unit/geom/CompoundCurveTest.cpp @@ -85,7 +85,8 @@ void object::test<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? - ensure(!cc_->getEnvelopeInternal()->isNull()); + geos::geom::Envelope expected(0, 2, 0, 2); + ensure("getEnvelopeInternal", cc_->getEnvelopeInternal()->equals(&expected)); // Geometry dimension functions ensure_equals("getDimension", cc_->getDimension(), geos::geom::Dimension::L); diff --git a/tests/unit/geom/CurvePolygonTest.cpp b/tests/unit/geom/CurvePolygonTest.cpp index 135f2dec35..0025235d9d 100644 --- a/tests/unit/geom/CurvePolygonTest.cpp +++ b/tests/unit/geom/CurvePolygonTest.cpp @@ -95,7 +95,15 @@ void object::test<2>() ensure_equals("getNumGeometries", cp_->getNumGeometries(), 1u); ensure_equals("getNumPoints", cp_->getNumPoints(), 14u); ensure_equals("getNumInteriorRing", cp_->getNumInteriorRing(), 1u); - ensure(!cp_->getEnvelopeInternal()->isNull()); + { + geos::geom::Envelope expected(0, 4, -0.618033988749895, 5); + const geos::geom::Envelope& actual = *cp_->getEnvelopeInternal(); + + ensure_equals("getEnvelopeInternal MinX", actual.getMinX(), expected.getMinX()); + ensure_equals("getEnvelopeInternal MinY", actual.getMinY(), expected.getMinY()); + ensure_equals("getEnvelopeInternal MaxX", actual.getMaxX(), expected.getMaxX()); + ensure_equals("getEnvelopeInternal MaxY", actual.getMaxY(), expected.getMaxY()); + } // Geometry dimension functions ensure_equals("getDimension", cp_->getDimension(), geos::geom::Dimension::A);