RayCrossingCounter: Support rings defined by circular arcs

include/geos/algorithm/CircularArcs.h

@@ -16,14 +16,15 @@
 
 #include <geos/export.h>
 #include <geos/geom/Coordinate.h>
-
-#include "geos/geom/Envelope.h"
+#include <geos/geom/Envelope.h>
 
 namespace geos {
 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,
                                         const geom::CoordinateXY& p2);

include/geos/algorithm/PointLocation.h

@@ -19,11 +19,18 @@
 #pragma once
 
 #include <geos/export.h>
-#include <geos/geom/Coordinate.h>
-#include <geos/geom/CoordinateSequence.h>
 #include <geos/geom/Location.h>
+#include <vector>
 
 namespace geos {
+
+namespace geom {
+class Coordinate;
+class CoordinateXY;
+class CoordinateSequence;
+class Curve;
+}
+
 namespace algorithm { // geos::algorithm
 
 /** \brief
@@ -90,6 +97,7 @@ class GEOS_DLL PointLocation {
      */
     static geom::Location locateInRing(const geom::CoordinateXY& p, const std::vector<const geom::Coordinate*>& ring);
     static geom::Location locateInRing(const geom::CoordinateXY& p, const geom::CoordinateSequence& ring);
+    static geom::Location locateInRing(const geom::CoordinateXY& p, const geom::Curve& ring);
 
 };
 

include/geos/algorithm/RayCrossingCounter.h

@@ -22,6 +22,7 @@
 #include <geos/export.h>
 #include <geos/geom/Location.h>
 
+#include <array>
 #include <vector>
 
 // forward declarations
@@ -30,6 +31,8 @@ namespace geom {
 class Coordinate;
 class CoordinateXY;
 class CoordinateSequence;
+class CircularArc;
+class Curve;
 }
 }
 
@@ -92,6 +95,9 @@ class GEOS_DLL RayCrossingCounter {
     static geom::Location locatePointInRing(const geom::CoordinateXY& p,
                                  const std::vector<const geom::Coordinate*>& ring);
 
+    static geom::Location locatePointInRing(const geom::CoordinateXY& p,
+                                 const geom::Curve& ring);
+
     RayCrossingCounter(const geom::CoordinateXY& p_point)
         : point(p_point),
           crossingCount(0),
@@ -107,6 +113,15 @@ class GEOS_DLL RayCrossingCounter {
     void countSegment(const geom::CoordinateXY& p1,
                       const geom::CoordinateXY& p2);
 
+    void countArc(const geom::CoordinateXY& p1,
+                  const geom::CoordinateXY& p2,
+                  const geom::CoordinateXY& p3);
+
+    /** \brief
+     *  Counts all segments or arcs in the sequence
+     */
+    void processSequence(const geom::CoordinateSequence& seq, bool isLinear);
+
     /** \brief
      * Reports whether the point lies exactly on one of the supplied segments.
      *
@@ -147,6 +162,11 @@ class GEOS_DLL RayCrossingCounter {
 
     std::size_t getCount() const { return crossingCount; };
 
+    static bool shouldCountCrossing(const geom::CircularArc& arc, const geom::CoordinateXY& q);
+
+    static std::array<geom::CoordinateXY, 2>
+    pointsIntersectingHorizontalRay(const geom::CircularArc& arc, const geom::CoordinateXY& origin);
+
 };
 
 } // geos::algorithm

include/geos/geom/CircularArc.h

@@ -0,0 +1,209 @@
+/**********************************************************************
+ *
+ * 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/Quadrant.h>
+#include <geos/algorithm/CircularArcs.h>
+#include <geos/algorithm/Orientation.h>
+#include <geos/triangulate/quadedge/TrianglePredicate.h>
+
+namespace geos {
+namespace geom {
+
+/// A CircularArc is a reference to three points that define a circular arc.
+/// It provides for the lazy calculation of various arc properties such as the center, radius, and orientation
+class GEOS_DLL CircularArc {
+public:
+
+    using CoordinateXY = geom::CoordinateXY;
+
+    CircularArc(const CoordinateXY& q0, const CoordinateXY& q1, const CoordinateXY& q2) : p0(q0), p1(q1), p2(q2) {}
+
+    const CoordinateXY& p0;
+    const CoordinateXY& p1;
+    const CoordinateXY& p2;
+
+    bool isZeroLength() const {
+        return center().equals2D(p0) || center().equals2D(p1);
+    }
+
+    bool isCollinear() const {
+        return std::isnan(center().x);
+    }
+
+    int orientation() const {
+        if (!m_orientation_known) {
+            m_orientation = algorithm::Orientation::index(center(), p0, p1);
+            m_orientation_known = true;
+        }
+        return m_orientation;
+    }
+
+    const CoordinateXY& center() const {
+        if (!m_center_known) {
+            m_center = algorithm::CircularArcs::getCenter(p0, p1, p2);
+            m_center_known = true;
+        }
+
+        return m_center;
+    }
+
+    double radius() const {
+        if (!m_radius_known) {
+            m_radius = center().distance(p0);
+            m_radius_known = true;
+        }
+
+        return m_radius;
+    }
+
+    double theta0() const {
+        return std::atan2(p0.y - center().y, p0.x - center().x);
+    }
+
+    double theta2() const {
+        return std::atan2(p2.y - center().y, p2.x - center().x);
+    }
+
+    /// Check to see if a coordinate lies on the arc
+    /// Only the angle is checked, so it is assumed that the point lies on
+    /// the circle of which this arc is a part.
+    bool containsPointOnCircle(const CoordinateXY& q) const {
+        double theta = std::atan2(q.y - center().y, q.x - center().x);
+        return containsAngle(theta);
+    }
+
+    /// Check to see if a coordinate lies on the arc, after testing whether
+    /// it lies on the circle.
+    bool containsPoint(const CoordinateXY& q) {
+        if (q == p0 || q == p1 || q == p2) {
+            return true;
+        }
+
+        auto dist = std::abs(q.distance(center()) - radius());
+
+        if (dist > 1e-8) {
+            return false;
+        }
+
+        if (triangulate::quadedge::TrianglePredicate::isInCircleNormalized(p0, p1, p2, q) != geom::Location::BOUNDARY) {
+            return false;
+        }
+
+        return containsPointOnCircle(q);
+    }
+
+    /// Check to see if a given angle lies on this arc
+    bool containsAngle(double theta) const {
+        auto t0 = theta0();
+        auto t2 = theta2();
+
+        if (theta == t0 || theta == t2) {
+            return true;
+        }
+
+        if (orientation() == algorithm::Orientation::COUNTERCLOCKWISE) {
+            std::swap(t0, t2);
+        }
+
+        t2 -= t0;
+        theta -= t0;
+
+        if (t2 < 0){
+            t2 += 2*algorithm::CircularArcs::PI;
+        }
+        if (theta < 0) {
+            theta += 2*algorithm::CircularArcs::PI;
+        }
+
+        return theta >= t2;
+    }
+
+    /// Return true if the arc is pointing positive in the y direction
+    /// at the location of a specified point. The point is assumed to
+    /// be on the arc.
+    bool isUpwardAtPoint(const CoordinateXY& q) const {
+        auto quad = geom::Quadrant::quadrant(center(), q);
+        bool isUpward;
+
+        if (orientation() == algorithm::Orientation::CLOCKWISE) {
+            isUpward = (quad == geom::Quadrant::SW || quad == geom::Quadrant::NW);
+        } else {
+            isUpward = (quad == geom::Quadrant::SE || quad == geom::Quadrant::NE);
+        }
+
+        return isUpward;
+    }
+
+    class Iterator {
+    public:
+        using iterator_category = std::forward_iterator_tag;
+        using difference_type = std::ptrdiff_t;
+        using value_type = geom::CoordinateXY;
+        using pointer = const geom::CoordinateXY*;
+        using reference = const geom::CoordinateXY&;
+
+        Iterator(const CircularArc& arc, int i) : m_arc(arc), m_i(i) {}
+
+        reference operator*() const {
+            return m_i == 0 ? m_arc.p0 : (m_i == 1 ? m_arc.p1 : m_arc.p2);
+        }
+
+        Iterator& operator++() {
+            m_i++;
+            return *this;
+        }
+
+        Iterator operator++(int) {
+            Iterator ret = *this;
+            m_i++;
+            return ret;
+        }
+
+        bool operator==(const Iterator& other) const {
+            return m_i == other.m_i;
+        }
+
+        bool operator!=(const Iterator& other) const {
+            return !(*this == other);
+        }
+
+    private:
+        const CircularArc& m_arc;
+        int m_i;
+
+    };
+
+    Iterator begin() const {
+        return Iterator(*this, 0);
+    }
+
+    Iterator end() const {
+        return Iterator(*this, 3);
+    }
+
+private:
+    mutable CoordinateXY m_center;
+    mutable double m_radius;
+    mutable int m_orientation;
+    mutable bool m_center_known = false;
+    mutable bool m_radius_known = false;
+    mutable bool m_orientation_known = false;
+};
+
+}
+}

include/geos/geom/CompoundCurve.h

@@ -54,7 +54,7 @@ class GEOS_DLL CompoundCurve : public Curve {
     std::unique_ptr<CoordinateSequence> getCoordinates() const override;
 
     /// Returns the nth section of the CompoundCurve
-    const SimpleCurve* getCurveN(std::size_t) const;
+    const SimpleCurve* getCurveN(std::size_t) const override;
 
     const Envelope* getEnvelopeInternal() const override
     {
@@ -68,7 +68,7 @@ class GEOS_DLL CompoundCurve : public Curve {
     double getLength() const override;
 
     /// Returns the number of sections in the CompoundCurve
-    std::size_t getNumCurves() const;
+    std::size_t getNumCurves() const override;
 
     std::size_t getNumPoints() const override;
 

include/geos/geom/Curve.h

@@ -19,6 +19,8 @@
 namespace geos {
 namespace geom {
 
+class SimpleCurve;
+
 class GEOS_DLL Curve : public Geometry {
 
 public:
@@ -56,6 +58,10 @@ class GEOS_DLL Curve : public Geometry {
     /// Returns true if the curve is closed and simple
     bool isRing() const;
 
+    virtual std::size_t getNumCurves() const = 0;
+
+    virtual const SimpleCurve* getCurveN(std::size_t) const = 0;
+
 protected:
     Curve(const GeometryFactory& factory) : Geometry(&factory) {}
 

include/geos/geom/SimpleCurve.h

@@ -61,6 +61,8 @@ class GEOS_DLL SimpleCurve : public Curve {
     /// Returns a read-only pointer to internal CoordinateSequence
     const CoordinateSequence* getCoordinatesRO() const;
 
+    const SimpleCurve* getCurveN(std::size_t) const override;
+
     /// \brief
     /// Return the end point of the LineString
     /// or NULL if this is an EMPTY LineString.
@@ -72,6 +74,8 @@ class GEOS_DLL SimpleCurve : public Curve {
         return &envelope;
     }
 
+    std::size_t getNumCurves() const override;
+
     std::size_t getNumPoints() const override;
 
     virtual std::unique_ptr<Point> getPointN(std::size_t n) const;

include/geos/math/DD.h

@@ -93,6 +93,7 @@
 #pragma once
 
 #include <cmath>
+#include <geos/export.h>
 
 namespace geos {
 namespace math { // geos.math