Add MaximumInscribedCircle fast exact calculation for simple shapes, port of locationtech/jts#1123

Author: pramsey

include/geos/algorithm/Angle.h

@@ -147,6 +147,19 @@ class GEOS_DLL Angle {
                                        const geom::CoordinateXY& tail,
                                        const geom::CoordinateXY& tip2);
 
+    /// Computes the angle of the unoriented bisector
+    /// of the smallest angle between two vectors.
+    /// The computed angle will be in the range (-Pi, Pi].
+    ///
+    /// @param tip1 the tip of v1
+    /// @param tail the tail of each vector
+    /// @param tip2 the tip of v2
+    /// @return the angle of the bisector between v1 and v2
+    ///
+    static double bisector(const geom::CoordinateXY& tip1,
+                           const geom::CoordinateXY& tail,
+                           const geom::CoordinateXY& tip2);
+
     /// Computes the interior angle between two segments of a ring.
     ///
     /// The ring is assumed to be oriented in a clockwise direction.
@@ -230,14 +243,23 @@ class GEOS_DLL Angle {
     /// @param rCos the result of cos(ang)
     ///
     static inline void sinCosSnap(const double ang, double& rSin, double& rCos) {
-        // calculate both; may be optimized with FSINCOS instruction
         rSin = std::sin(ang);
         rCos = std::cos(ang);
         // snap near-zero values
         if (std::fabs(rSin) < 5e-16) rSin = 0.0;
         if (std::fabs(rCos) < 5e-16) rCos = 0.0;
     }
 
+    /// Projects a point by a given angle and distance.
+    ///
+    /// @param p the point to project
+    /// @param angle the angle at which to project
+    /// @param dist the distance to project
+    /// @return the projected point
+    ///
+    static geom::CoordinateXY project(const geom::CoordinateXY& p,
+        double angle, double dist);
+
 };
 
 

include/geos/algorithm/construct/ExactMaxInscribedCircle.h

@@ -0,0 +1,122 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (c) 2020 Martin Davis
+ * Copyright (C) 2020 Paul Ramsey <pramsey@cleverelephant.ca>
+ *
+ * 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.
+ *
+ **********************************************************************
+ *
+ * Last port: algorithm/construct/ExactMaxInscribedCircle.java
+ * https://github.com/locationtech/jts/commit/8d5ced1b784d232e1eecf5df4b71873e8822336a
+ *
+ **********************************************************************/
+
+#pragma once
+
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/LineSegment.h>
+
+#include <array>
+// #include <queue>
+
+
+
+namespace geos {
+namespace geom {
+class CoordinateSequence;
+class Geometry;
+class Polygon;
+}
+}
+
+namespace geos {
+namespace algorithm { // geos::algorithm
+namespace construct { // geos::algorithm::construct
+
+/**
+ * Computes the Maximum Inscribed Circle for some kinds of convex polygons.
+ * It determines the circle center point by computing Voronoi node points
+ * and testing them for distance to generating edges.
+ * This is more precise than iterated approximation,
+ * and faster for small polygons (such as triangles and convex quadrilaterals).
+ *
+ * @author Martin Davis
+ *
+ */
+class GEOS_DLL ExactMaxInscribedCircle {
+
+    using CoordinateSequence = geos::geom::CoordinateSequence;
+    using CoordinateXY = geos::geom::CoordinateXY;
+    using LineSegment = geos::geom::LineSegment;
+    using Polygon = geos::geom::Polygon;
+
+public:
+
+    ExactMaxInscribedCircle();
+
+    /**
+     * Tests whether a given geometry is supported by this class.
+     * Currently only triangles and convex quadrilaterals are supported.
+     *
+     * @param geom an areal geometry
+     * @return true if the geometry shape can be evaluated
+     */
+    static bool isSupported(const geom::Geometry* geom);
+
+    static std::pair<CoordinateXY, CoordinateXY> computeRadius(const Polygon* polygon);
+
+
+private:
+
+    static bool isTriangle(const Polygon* polygon);
+
+    static bool isQuadrilateral(const Polygon* polygon);
+
+    static std::pair<CoordinateXY, CoordinateXY> computeTriangle(const CoordinateSequence* ring);
+
+    /**
+     * The Voronoi nodes of a convex polygon occur at the intersection point
+     * of two bisectors of each triplet of edges.
+     * The Maximum Inscribed Circle center is the node
+     * is the farthest distance from the generating edges.
+     * For a quadrilateral there are 4 distinct edge triplets,
+     * at each edge with its adjacent edges.
+     *
+     * @param ring the polygon ring
+     * @return an array containing the incircle center and radius points
+     */
+    static std::pair<CoordinateXY, CoordinateXY> computeConvexQuadrilateral(const CoordinateSequence* ring);
+
+    static std::array<LineSegment, 4> computeBisectors(const CoordinateSequence* ptsCW,
+                                                       double diameter);
+
+    static CoordinateXY nearestEdgePt(const CoordinateSequence* ring,
+                                      const CoordinateXY& pt);
+
+    static LineSegment computeConvexBisector(const CoordinateSequence* pts,
+                                             std::size_t index, double len);
+
+    static bool isConvex(const Polygon* polygon);
+
+    static bool isConvex(const CoordinateSequence* ring);
+
+    static bool isConvex(const CoordinateXY& p0,
+                         const CoordinateXY& p1,
+                         const CoordinateXY& p2);
+
+    static bool isPointInConvexRing(const CoordinateSequence* ringCW,
+                                const CoordinateXY& p);
+
+};
+
+
+} // geos::algorithm::construct
+} // geos::algorithm
+} // geos

include/geos/algorithm/construct/MaximumInscribedCircle.h

@@ -139,6 +139,8 @@ class GEOS_DLL MaximumInscribedCircle {
     double distanceToBoundary(const geom::Point& pt);
     double distanceToBoundary(double x, double y);
     void compute();
+    void computeApproximation();
+    void createResult(const geom::CoordinateXY& c, const geom::CoordinateXY& r);
 
     /* private class */
     class Cell {

include/geos/geom/Triangle.h

@@ -45,6 +45,11 @@ class GEOS_DLL Triangle {
      */
     void inCentre(CoordinateXY& resultPoint);
 
+    static CoordinateXY inCentre(
+        const CoordinateXY& p0,
+        const CoordinateXY& p1,
+        const CoordinateXY& p2);
+
     /** \brief
      * Computes the circumcentre of a triangle.
      *

src/algorithm/Angle.cpp

@@ -118,6 +118,19 @@ Angle::angleBetweenOriented(const geom::CoordinateXY& tip1,
     return angDel;
 }
 
+
+/* public static */
+double
+Angle::bisector(const geom::CoordinateXY& tip1,
+                const geom::CoordinateXY& tail,
+                const geom::CoordinateXY& tip2)
+{
+    double angDel = angleBetweenOriented(tip1, tail, tip2);
+    double angBi = angle(tail, tip1) + angDel / 2;
+    return normalize(angBi);
+}
+
+
 /* public static */
 double
 Angle::interiorAngle(const geom::CoordinateXY& p0, const geom::CoordinateXY& p1,
@@ -187,20 +200,32 @@ Angle::diff(double ang1, double ang2)
 {
     double delAngle;
 
-    if(ang1 < ang2) {
+    if (ang1 < ang2) {
         delAngle = ang2 - ang1;
     }
     else {
         delAngle = ang1 - ang2;
     }
 
-    if(delAngle > MATH_PI) {
+    if (delAngle > MATH_PI) {
         delAngle = PI_TIMES_2 - delAngle;
     }
 
     return delAngle;
 }
 
+/* public static */
+geom::CoordinateXY
+Angle::project(const geom::CoordinateXY& p, double angle, double dist)
+{
+    double cosSnap, sinSnap;
+    sinCosSnap(angle, sinSnap, cosSnap);
+    double x = p.x + dist * cosSnap;
+    double y = p.y + dist * sinSnap;
+    return geom::CoordinateXY(x, y);
+}
+
+
 } // namespace geos.algorithm
 } //namespace geos