WIP: Add CircularArcIntersector

include/geos/algorithm/Angle.h

@@ -204,6 +204,15 @@ class GEOS_DLL Angle {
     ///
     static double normalizePositive(double angle);
 
+    /// Returns true if angle x is within the counterclockwise
+    /// arc from angle a to angle b
+    ///
+    /// @param angle angle to test
+    /// @param from starting angle of arc
+    /// @param to ending angle of arc
+    ///
+    /// @return true if `angle` is within [from, to]
+    static bool isWithinCCW(double angle, double from, double to);
 
     /// Computes the unoriented smallest difference between two angles.
     ///

include/geos/algorithm/CircularArcIntersector.h

@@ -0,0 +1,97 @@
+/**********************************************************************
+ *
+ * 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 <array>
+#include <cstdint>
+
+#include <geos/export.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/CircularArc.h>
+#include <geos/geom/LineSegment.h>
+
+namespace geos::algorithm {
+
+class GEOS_DLL CircularArcIntersector {
+public:
+    using CoordinateXY = geom::CoordinateXY;
+    using CircularArc = geom::CircularArc;
+    using Envelope = geom::Envelope;
+
+    enum intersection_type : uint8_t {
+        NO_INTERSECTION = 0,
+        ONE_POINT_INTERSECTION = 1,
+        TWO_POINT_INTERSECTION = 2,
+        COCIRCULAR_INTERSECTION = 3,
+    };
+
+    intersection_type getResult() const
+    {
+        return result;
+    }
+
+    const CoordinateXY& getPoint(std::uint8_t i) const
+    {
+        return intPt[i];
+    }
+
+    const std::array<CoordinateXY, 3>& getArc(std::uint8_t i) const
+    {
+        return intArc[i];
+    }
+
+    std::uint8_t getNumPoints() const
+    {
+        return nPt;
+    }
+
+    std::uint8_t getNumArcs() const
+    {
+        return nArc;
+    }
+
+    /// Determines whether and where a circular arc intersects a line segment.
+    ///
+    /// Sets the appropriate value of intersection_type and stores the intersection
+    /// points, if any.
+    void intersects(const CircularArc& arc, const CoordinateXY& p0, const CoordinateXY& p1);
+
+    void intersects(const CircularArc& arc, const geom::LineSegment& seg)
+    {
+        intersects(arc, seg.p0, seg.p1);
+    }
+
+    /// Determines whether and where two circular arcs intersect.
+    ///
+    ///	Sets the appropriate value of intesection_type and stores the intersection
+    /// points and/or arcs, if any.
+    void intersects(const CircularArc& arc1, const CircularArc& arc2);
+
+    static int
+    circleIntersects(const CoordinateXY& center, double r, const CoordinateXY& p0, const CoordinateXY& p1, CoordinateXY& isect0, CoordinateXY& isect1);
+
+private:
+
+    void intersects(const CoordinateXY& p0, const CoordinateXY& p1, const CoordinateXY& q0, const CoordinateXY& q1);
+
+    std::array<CoordinateXY, 2> intPt;
+    std::array<std::array<CoordinateXY, 3>, 2> intArc;
+    intersection_type result;
+    std::uint8_t nPt = 0;
+    std::uint8_t nArc = 0;
+
+};
+
+}

include/geos/algorithm/CircularArcs.h

@@ -14,6 +14,8 @@
 
 #pragma once
 
+#include <array>
+
 #include <geos/export.h>
 #include <geos/geom/Coordinate.h>
 #include <geos/geom/Envelope.h>
@@ -31,6 +33,15 @@ class GEOS_DLL CircularArcs {
     /// 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);
+
+
+    /// Return three points defining a arc defined by a circle center, radius, and start/end angles
+    static std::array<geom::CoordinateXY, 3>
+    createArc(const geom::CoordinateXY& center, double radius, double start, double end, bool ccw);
+
+    /// Return the point defined by a circle center, radius, and angle
+    static geom::CoordinateXY createPoint(const geom::CoordinateXY& center, double radius, double theta);
+
 };
 
 }

include/geos/geom/CircularArc.h

@@ -83,7 +83,7 @@ class GEOS_DLL CircularArc {
 
     /// Returns whether this arc forms a straight line (p0, p1, and p2 are collinear)
     bool isLinear() const {
-        return std::isnan(getRadius());
+        return !std::isfinite(getRadius());
     }
 
     /// Return the inner angle of the sector associated with this arc

include/geos/math/DD.h

@@ -158,6 +158,7 @@ class GEOS_DLL DD {
         friend GEOS_DLL DD operator* (const DD &lhs, double rhs);
         friend GEOS_DLL DD operator/ (const DD &lhs, const DD &rhs);
         friend GEOS_DLL DD operator/ (const DD &lhs, double rhs);
+        friend GEOS_DLL DD operator- (const DD& x);
 
         static DD determinant(const DD &x1, const DD &y1, const DD &x2, const DD &y2);
         static DD determinant(double x1, double y1, double x2, double y2);

src/algorithm/Angle.cpp

@@ -181,6 +181,15 @@ Angle::normalizePositive(double angle)
     return angle;
 }
 
+bool
+Angle::isWithinCCW(double x, double a, double b) {
+    if (b > a) {
+        return x >= a && x <= b;
+    } else {
+        return x >= a || x <= b;
+    }
+}
+
 /* public static */
 double
 Angle::diff(double ang1, double ang2)