Implement getEnvelopeInternal for curved types

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 <geos/export.h>
+#include <geos/geom/Coordinate.h>
+
+#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);
+};
+
+}
+}

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<CoordinateSequence> && pts,
+                   const GeometryFactory& newFactory);
+
+    void geometryChangedAction() override
+    {
+        envelope = computeEnvelopeInternal(false);
+    }
+
 };
 
 

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();

include/geos/geom/SimpleCurve.h

@@ -126,14 +126,10 @@ class GEOS_DLL SimpleCurve : public Curve {
     SimpleCurve(const SimpleCurve& other);
 
     SimpleCurve(std::unique_ptr<CoordinateSequence>&& 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<CoordinateSequence> points;

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 <geos/algorithm/CircularArcs.h>
+
+#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;
+        }
+    }
+}
+
+}
+}

src/geom/CircularString.cpp

@@ -20,6 +20,15 @@
 namespace geos {
 namespace geom {
 
+/*public*/
+CircularString::CircularString(std::unique_ptr<CoordinateSequence> && newCoords,
+                               const GeometryFactory& factory)
+    :
+    SimpleCurve(std::move(newCoords), false, factory)
+{
+    //validateConstruction();
+}
+
 CircularString::~CircularString() = default;
 
 std::string

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();
 }

src/geom/SimpleCurve.cpp

@@ -17,6 +17,7 @@
 
 #include <geos/geom/SimpleCurve.h>
 
+#include <geos/algorithm/CircularArcs.h>
 #include <geos/algorithm/Orientation.h>
 #include <geos/geom/CoordinateFilter.h>
 #include <geos/geom/GeometryFactory.h>
@@ -38,21 +39,33 @@ SimpleCurve::SimpleCurve(const SimpleCurve& other)
 }
 
 SimpleCurve::SimpleCurve(std::unique_ptr<CoordinateSequence>&& newCoords,
+                         bool isLinear,
             const GeometryFactory& factory)
     : Curve(factory),
     points(newCoords ? std::move(newCoords) : std::make_unique<CoordinateSequence>()),
-    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<CoordinateXY>(i-2),
+                                                    points->getAt<CoordinateXY>(i-1),
+                                                    points->getAt<CoordinateXY>(i));
+        }
+        return e;
+    }
 }
 
 std::unique_ptr<CoordinateSequence>