CircularArc.h

/**********************************************************************
 *
 * 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/LineSegment.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 {
 
class GEOS_DLL CircularArc {
public:
 
    using CoordinateXY = geom::CoordinateXY;
 
    CircularArc() : CircularArc({0, 0}, {0, 0}, {0, 0}) {}
 
    CircularArc(const CoordinateXY& q0, const CoordinateXY& q1, const CoordinateXY& q2)
        : p0(q0)
        , p1(q1)
        , p2(q2)
        , m_center_known(false)
        , m_radius_known(false)
        , m_orientation_known(false)
    {}
 
    CircularArc(double theta0, double theta2, const CoordinateXY& center, double radius, int orientation)
        : p0(algorithm::CircularArcs::createPoint(center, radius, theta0)),
          p1(algorithm::CircularArcs::createPoint(center, radius, algorithm::CircularArcs::getMidpointAngle(theta0, theta2, orientation==algorithm::Orientation::COUNTERCLOCKWISE))),
          p2(algorithm::CircularArcs::createPoint(center, radius, theta2)),
          m_center(center),
          m_radius(radius),
          m_orientation(orientation),
          m_center_known(true),
          m_radius_known(true),
          m_orientation_known(true)
    {}
 
    CircularArc(const CoordinateXY& q0, const CoordinateXY& q2, const CoordinateXY& center, double radius, int orientation)
        : p0(q0),
          p1(algorithm::CircularArcs::getMidpoint(q0, q2, center, radius, orientation==algorithm::Orientation::COUNTERCLOCKWISE)),
          p2(q2),
          m_center(center),
          m_radius(radius),
          m_orientation(orientation),
          m_center_known(true),
          m_radius_known(true),
          m_orientation_known(true)
    {}
 
    CoordinateXY p0;
    CoordinateXY p1;
    CoordinateXY p2;
 
    int getOrientation() const {
        if (!m_orientation_known) {
            m_orientation = algorithm::Orientation::index(p0, p1, p2);
            m_orientation_known = true;
        }
        return m_orientation;
    }
 
    bool isCCW() const {
        return getOrientation() == algorithm::Orientation::COUNTERCLOCKWISE;
    }
 
    const CoordinateXY& getCenter() const {
        if (!m_center_known) {
            m_center = algorithm::CircularArcs::getCenter(p0, p1, p2);
            m_center_known = true;
        }
 
        return m_center;
    }
 
    double getRadius() const {
        if (!m_radius_known) {
            m_radius = getCenter().distance(p0);
            m_radius_known = true;
        }
 
        return m_radius;
    }
 
    bool isCircle() const {
        return p0.equals(p2);
    }
 
    bool isLinear() const {
        return !std::isfinite(getRadius());
    }
 
    double getAngle() const {
        if (isCircle()) {
            return 2*MATH_PI;
        }
 
        auto t0 = theta0();
        auto t2 = theta2();
 
        if (getOrientation() == algorithm::Orientation::COUNTERCLOCKWISE) {
            std::swap(t0, t2);
        }
 
        if (t0 < t2) {
            t0 += 2*MATH_PI;
        }
 
        auto diff = t0-t2;
 
        return diff;
    }
 
    double getLength() const {
        if (isLinear()) {
            return p0.distance(p2);
        }
 
        return getAngle()*getRadius();
    }
 
    double getArea() const {
        if (isLinear()) {
            return 0;
        }
 
        auto R = getRadius();
        auto theta = getAngle();
        return R*R/2*(theta - std::sin(theta));
    }
 
    double getSagitta() const {
        CoordinateXY midpoint = algorithm::CircularArcs::getMidpoint(p0, p2, getCenter(), getRadius(), isCCW());
        return algorithm::Distance::pointToSegment(midpoint, p0, p2);
    }
 
    double theta0() const {
        return algorithm::CircularArcs::getAngle(p0, getCenter());
    }
 
    double theta2() const {
        return algorithm::CircularArcs::getAngle(p2, getCenter());
    }
 
    bool containsPointOnCircle(const CoordinateXY& q) const {
        double theta = std::atan2(q.y - getCenter().y, q.x - getCenter().x);
        return containsAngle(theta);
    }
 
    bool containsPoint(const CoordinateXY& q) const {
        if (q == p0 || q == p1 || q == p2) {
            return true;
        }
 
        //auto dist = std::abs(q.distance(getCenter()) - getRadius());
 
        //if (dist > 1e-8) {
        //    return false;
        //}
 
        if (triangulate::quadedge::TrianglePredicate::isInCircleRobust(p0, p1, p2, q) != geom::Location::BOUNDARY) {
            return false;
        }
 
        return containsPointOnCircle(q);
    }
 
    bool containsAngle(double theta) const {
        auto t0 = theta0();
        auto t2 = theta2();
 
        if (theta == t0 || theta == t2) {
            return true;
        }
 
        if (getOrientation() == algorithm::Orientation::COUNTERCLOCKWISE) {
            std::swap(t0, t2);
        }
 
        t2 -= t0;
        theta -= t0;
 
        if (t2 < 0) {
            t2 += 2*MATH_PI;
        }
        if (theta < 0) {
            theta += 2*MATH_PI;
        }
 
        return theta >= t2;
    }
 
    bool isUpwardAtPoint(const CoordinateXY& q) const {
        auto quad = geom::Quadrant::quadrant(getCenter(), q);
        bool isUpward;
 
        if (getOrientation() == algorithm::Orientation::CLOCKWISE) {
            isUpward = (quad == geom::Quadrant::SW || quad == geom::Quadrant::NW);
        } else {
            isUpward = (quad == geom::Quadrant::SE || quad == geom::Quadrant::NE);
        }
 
        return isUpward;
    }
 
    // Split an arc at a specified point.
    // The point is assumed to be o the arc.
    std::pair<CircularArc, CircularArc> splitAtPoint(const CoordinateXY& q) const {
        return {
            CircularArc(p0, q, getCenter(), getRadius(), getOrientation()),
            CircularArc(q, p2, getCenter(), getRadius(), getOrientation())
        };
    }
 
    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;
};
 
}
}