Add Angle::sinCosSnap to avoid small errors, e.g. with buffer operations

NEWS.md

@@ -2,6 +2,7 @@
 20xx-xx-xx
 
 - New things:
+  - Add Angle::SinCos to avoid small errors, e.g. with buffer operations (GH-978, Mike Taves)
 
 - Breaking Changes
 

include/geos/algorithm/Angle.h

@@ -216,6 +216,27 @@ class GEOS_DLL Angle {
     ///         (in range [0, Pi] )
     ///
     static double diff(double ang1, double ang2);
+
+    /// \brief
+    /// Computes both sin and cos of an angle.
+    ///
+    /// The angle does not need to be normalized. Unlike std::sin
+    /// and std::cos, this method will clip near-zero values to zero
+    /// for (e.g.) sin(pi) and cos(pi/2).
+    ///
+    /// @param ang the input angle (in radians)
+    /// @param rSin the result of sin(ang)
+    /// @param rCos the result of cos(ang)
+    ///
+    static inline void SinCos(const double ang, double& rSin, double& rCos) {
+        // calculate both; may be optimized with FSINCOS instruction
+        rSin = std::sin(ang);
+        rCos = std::cos(ang);
+        // clip near zero values
+        if (std::fabs(rSin) < 5e-16) rSin = 0.0;
+        if (std::fabs(rCos) < 5e-16) rCos = 0.0;
+    }
+
 };
 
 

src/algorithm/Angle.cpp

@@ -195,7 +195,7 @@ Angle::diff(double ang1, double ang2)
     }
 
     if(delAngle > MATH_PI) {
-        delAngle = (2 * MATH_PI) - delAngle;
+        delAngle = PI_TIMES_2 - delAngle;
     }
 
     return delAngle;

src/operation/buffer/BufferParameters.cpp

@@ -19,6 +19,7 @@
 #include <cstdlib> // for std::abs()
 #include <cmath> // for cos
 
+#include <geos/algorithm/Angle.h>
 #include <geos/constants.h>
 #include <geos/operation/buffer/BufferParameters.h>
 
@@ -95,7 +96,7 @@ BufferParameters::setQuadrantSegments(int quadSegs)
 double
 BufferParameters::bufferDistanceError(int quadSegs)
 {
-    double alpha = MATH_PI / 2.0 / quadSegs;
+    double alpha = algorithm::Angle::PI_OVER_2 / quadSegs;
     return 1 - cos(alpha / 2.0);
 }
 

src/operation/buffer/OffsetSegmentGenerator.cpp

@@ -79,11 +79,11 @@ OffsetSegmentGenerator::OffsetSegmentGenerator(
 {
     // compute intersections in full precision, to provide accuracy
     // the points are rounded as they are inserted into the curve line
-    filletAngleQuantum = MATH_PI / 2.0 / bufParams.getQuadrantSegments();
+    filletAngleQuantum = Angle::PI_OVER_2 / bufParams.getQuadrantSegments();
 
     int quadSegs = bufParams.getQuadrantSegments();
     if (quadSegs < 1) quadSegs = 1;
-    filletAngleQuantum = MATH_PI / 2.0 / quadSegs;
+    filletAngleQuantum = Angle::PI_OVER_2 / quadSegs;
 
     /*
      * Non-round joins cause issues with short closing segments,
@@ -208,7 +208,7 @@ OffsetSegmentGenerator::addLineEndCap(const Coordinate& p0, const Coordinate& p1
     case BufferParameters::CAP_ROUND:
         // add offset seg points with a fillet between them
         segList.addPt(offsetL.p1);
-        addDirectedFillet(p1, angle + MATH_PI / 2.0, angle - MATH_PI / 2.0,
+        addDirectedFillet(p1, angle + Angle::PI_OVER_2, angle - Angle::PI_OVER_2,
                   Orientation::CLOCKWISE, distance);
         segList.addPt(offsetR.p1);
         break;
@@ -221,8 +221,10 @@ OffsetSegmentGenerator::addLineEndCap(const Coordinate& p0, const Coordinate& p1
         // add a square defined by extensions of the offset
         // segment endpoints
         Coordinate squareCapSideOffset;
-        squareCapSideOffset.x = fabs(distance) * cos(angle);
-        squareCapSideOffset.y = fabs(distance) * sin(angle);
+        double sinangle, cosangle;
+        Angle::SinCos(angle, sinangle, cosangle);
+        squareCapSideOffset.x = fabs(distance) * cosangle;
+        squareCapSideOffset.y = fabs(distance) * sinangle;
 
         Coordinate squareCapLOffset(
             offsetL.p1.x + squareCapSideOffset.x,
@@ -250,12 +252,12 @@ OffsetSegmentGenerator::addDirectedFillet(const Coordinate& p, const Coordinate&
 
     if(direction == Orientation::CLOCKWISE) {
         if(startAngle <= endAngle) {
-            startAngle += 2.0 * MATH_PI;
+            startAngle += Angle::PI_TIMES_2;
         }
     }
     else {    // direction==COUNTERCLOCKWISE
         if(startAngle >= endAngle) {
-            startAngle -= 2.0 * MATH_PI;
+            startAngle -= Angle::PI_TIMES_2;
         }
     }
 
@@ -277,13 +279,13 @@ OffsetSegmentGenerator::addDirectedFillet(const Coordinate& p, double startAngle
     // no segments because angle is less than increment-nothing to do!
     if(nSegs < 1) return;
 
-    // double initAngle, currAngleInc;
     double angleInc = totalAngle / nSegs;
+    double sinangle, cosangle;
     Coordinate pt;
     for (int i = 0; i < nSegs; i++) {
-        double angle = startAngle + directionFactor * i * angleInc;
-        pt.x = p.x + radius * cos(angle);
-        pt.y = p.y + radius * sin(angle);
+        Angle::SinCos(startAngle + directionFactor * i * angleInc, sinangle, cosangle);
+        pt.x = p.x + radius * cosangle;
+        pt.y = p.y + radius * sinangle;
         segList.addPt(pt);
     }
 }
@@ -295,7 +297,7 @@ OffsetSegmentGenerator::createCircle(const Coordinate& p, double p_distance)
     // add start point
     Coordinate pt(p.x + p_distance, p.y);
     segList.addPt(pt);
-    addDirectedFillet(p, 0.0, 2.0 * MATH_PI, -1, p_distance);
+    addDirectedFillet(p, 0.0, Angle::PI_TIMES_2, -1, p_distance);
     segList.closeRing();
 }
 
@@ -531,7 +533,7 @@ OffsetSegmentGenerator::addLimitedMitreJoin(
     Coordinate bevelMidPt = project(cornerPt, p_mitreLimitDistance, dirBisectorOut);
 
     // slope angle of bevel segment
-    double dirBevel = Angle::normalize(dirBisectorOut + MATH_PI/2.0);
+    double dirBevel = Angle::normalize(dirBisectorOut + Angle::PI_OVER_2);
 
     // compute the candidate bevel segment by projecting both sides of the midpoint
     Coordinate bevel0 = project(bevelMidPt, p_distance, dirBevel);
@@ -580,8 +582,10 @@ OffsetSegmentGenerator::extend(const LineSegment& seg, double dist)
 Coordinate
 OffsetSegmentGenerator::project(const Coordinate& pt, double d, double dir)
 {
-    double x = pt.x + d * std::cos(dir);
-    double y = pt.y + d * std::sin(dir);
+    double sindir, cosdir;
+    Angle::SinCos(dir, sindir, cosdir);
+    double x = pt.x + d * cosdir;
+    double y = pt.y + d * sindir;
     return Coordinate(x, y);
 }
 

src/util/GeometricShapeFactory.cpp

@@ -17,6 +17,7 @@
  *
  **********************************************************************/
 
+#include <geos/algorithm/Angle.h>
 #include <geos/constants.h>
 #include <geos/util/GeometricShapeFactory.h>
 #include <geos/geom/Coordinate.h>
@@ -31,6 +32,7 @@
 #include <memory>
 
 
+using namespace geos::algorithm;
 using namespace geos::geom;
 
 namespace geos {
@@ -134,10 +136,11 @@ GeometricShapeFactory::createCircle()
 
     auto pts = detail::make_unique<CoordinateSequence>(nPts + 1);
     uint32_t iPt = 0;
+    double sinang, cosang;
     for(uint32_t i = 0; i < nPts; i++) {
-        double ang = i * (2 * 3.14159265358979 / nPts);
-        double x = xRadius * cos(ang) + centreX;
-        double y = yRadius * sin(ang) + centreY;
+        Angle::SinCos(i * Angle::PI_TIMES_2 / nPts, sinang, cosang);
+        double x = xRadius * cosang + centreX;
+        double y = yRadius * sinang + centreY;
         (*pts)[iPt++] = coord(x, y);
     }
     (*pts)[iPt++] = (*pts)[0];
@@ -158,17 +161,18 @@ GeometricShapeFactory::createArc(double startAng, double angExtent)
     env.reset();
 
     double angSize = angExtent;
-    if(angSize <= 0.0 || angSize > 2 * MATH_PI) {
-        angSize = 2 * MATH_PI;
+    if(angSize <= 0.0 || angSize > Angle::PI_TIMES_2) {
+        angSize = Angle::PI_TIMES_2;
     }
     double angInc = angSize / (nPts - 1);
 
     auto pts = detail::make_unique<CoordinateSequence>(nPts);
     uint32_t iPt = 0;
+    double sinang, cosang;
     for(uint32_t i = 0; i < nPts; i++) {
-        double ang = startAng + i * angInc;
-        double x = xRadius * cos(ang) + centreX;
-        double y = yRadius * sin(ang) + centreY;
+        Angle::SinCos(startAng + i * angInc, sinang, cosang);
+        double x = xRadius * cosang + centreX;
+        double y = yRadius * sinang + centreY;
         (*pts)[iPt++] = coord(x, y);
     }
     auto line = geomFact->createLineString(std::move(pts));
@@ -187,18 +191,19 @@ GeometricShapeFactory::createArcPolygon(double startAng, double angExtent)
     env.reset();
 
     double angSize = angExtent;
-    if(angSize <= 0.0 || angSize > 2 * MATH_PI) {
-        angSize = 2 * MATH_PI;
+    if(angSize <= 0.0 || angSize > Angle::PI_TIMES_2) {
+        angSize = Angle::PI_TIMES_2;
     }
     double angInc = angSize / (nPts - 1);
 
     auto pts = detail::make_unique<CoordinateSequence>(nPts + 2);
     uint32_t iPt = 0;
     (*pts)[iPt++] = coord(centreX, centreY);
+    double sinang, cosang;
     for(uint32_t i = 0; i < nPts; i++) {
-        double ang = startAng + i * angInc;
-        double x = xRadius * cos(ang) + centreX;
-        double y = yRadius * sin(ang) + centreY;
+        Angle::SinCos(startAng + i * angInc, sinang, cosang);
+        double x = xRadius * cosang + centreX;
+        double y = yRadius * sinang + centreY;
         (*pts)[iPt++] = coord(x, y);
     }
     (*pts)[iPt++] = coord(centreX, centreY);

tests/unit/algorithm/AngleTest.cpp

@@ -162,6 +162,45 @@ void object::test<5>
         TOL);
 }
 
+// testSinCos()
+template<>
+template<>
+void object::test<6>
+()
+{
+    double rSin, rCos;
+
+    // -720 to 720 degrees with 1 degree increments
+    for (int angdeg = -720; angdeg <= 720; angdeg++) {
+        const double ang = Angle::toRadians(angdeg);
+
+        Angle::SinCos(ang, rSin, rCos);
+
+        double cSin = std::sin(ang);
+        double cCos = std::cos(ang);
+        if ( (angdeg % 90) == 0 ) {
+            // not always the same for multiples of 90 degrees
+            ensure(std::to_string(angdeg), std::fabs(rSin - cSin) < 1e-15);
+            ensure(std::to_string(angdeg), std::fabs(rCos - cCos) < 1e-15);
+        } else {
+            ensure_equals(std::to_string(angdeg), rSin, cSin);
+            ensure_equals(std::to_string(angdeg), rCos, cCos);
+        }
+
+    }
+
+    // use radian increments that don't snap to exact degrees or zero
+    for (double angrad = -6.3; angrad < 6.3; angrad += 0.013) {
+
+        Angle::SinCos(angrad, rSin, rCos);
+
+        ensure_equals(std::to_string(angrad), rSin, std::sin(angrad));
+        ensure_equals(std::to_string(angrad), rCos, std::cos(angrad));
+
+    }
+
+}
+
 
 } // namespace tut
 

tests/unit/capi/GEOSBufferTest.cpp

@@ -258,7 +258,7 @@ void object::test<9>
     ensure_area(area_, 150.0, 0.001);
 
     ensure_equals(std::string(wkt_), std::string(
-                      "POLYGON ((10.0000000000000000 15.0000000000000000, 15.0000000000000000 15.0000000000000000, 15.0000000000000000 5.0000000000000000, 5.0000000000000000 5.0000000000000000, 0.0000000000000000 5.0000000000000009, 0.0000000000000000 15.0000000000000000, 10.0000000000000000 15.0000000000000000))"
+                      "POLYGON ((10.0000000000000000 15.0000000000000000, 15.0000000000000000 15.0000000000000000, 15.0000000000000000 5.0000000000000000, 5.0000000000000000 5.0000000000000000, 0.0000000000000000 5.0000000000000000, 0.0000000000000000 15.0000000000000000, 10.0000000000000000 15.0000000000000000))"
                   ));
 }
 
@@ -311,7 +311,7 @@ void object::test<11>
     ensure_area(area_, 250.0, 0.001);
 
     ensure_equals(std::string(wkt_), std::string(
-                      "POLYGON ((5.0000000000000000 15.0000000000000000, 5.0000000000000000 20.0000000000000000, 5.0000000000000000 25.0000000000000000, 15.0000000000000000 25.0000000000000000, 15.0000000000000000 5.0000000000000000, 5.0000000000000000 5.0000000000000000, 0.0000000000000000 5.0000000000000009, 0.0000000000000000 15.0000000000000000, 5.0000000000000000 15.0000000000000000))"
+                      "POLYGON ((5.0000000000000000 15.0000000000000000, 5.0000000000000000 20.0000000000000000, 5.0000000000000000 25.0000000000000000, 15.0000000000000000 25.0000000000000000, 15.0000000000000000 5.0000000000000000, 5.0000000000000000 5.0000000000000000, 0.0000000000000000 5.0000000000000000, 0.0000000000000000 15.0000000000000000, 5.0000000000000000 15.0000000000000000))"
                   ));
 }
 
@@ -339,7 +339,7 @@ void object::test<12>
     ensure_area(area_, 237.5, 0.001);
 
     ensure_equals(std::string(wkt_), std::string(
-                      "POLYGON ((5.0000000000000000 15.0000000000000000, 5.0000000000000000 20.0000000000000000, 5.0000000000000000 25.0000000000000000, 15.0000000000000000 25.0000000000000000, 15.0000000000000000 10.0000000000000000, 10.0000000000000000 5.0000000000000000, 5.0000000000000000 5.0000000000000000, 0.0000000000000000 5.0000000000000009, 0.0000000000000000 15.0000000000000000, 5.0000000000000000 15.0000000000000000))"
+                      "POLYGON ((5.0000000000000000 15.0000000000000000, 5.0000000000000000 20.0000000000000000, 5.0000000000000000 25.0000000000000000, 15.0000000000000000 25.0000000000000000, 15.0000000000000000 10.0000000000000000, 10.0000000000000000 5.0000000000000000, 5.0000000000000000 5.0000000000000000, 0.0000000000000000 5.0000000000000000, 0.0000000000000000 15.0000000000000000, 5.0000000000000000 15.0000000000000000))"
                   ));
 }
 
@@ -368,7 +368,7 @@ void object::test<13>
     ensure_area(area_, 237.5, 0.001);
 
     ensure_equals(std::string(wkt_), std::string(
-                      "POLYGON ((5.0000000000000000 15.0000000000000000, 5.0000000000000000 20.0000000000000000, 5.0000000000000000 25.0000000000000000, 15.0000000000000000 25.0000000000000000, 15.0000000000000000 10.0000000000000000, 10.0000000000000000 5.0000000000000000, 5.0000000000000000 5.0000000000000000, 0.0000000000000000 5.0000000000000009, 0.0000000000000000 15.0000000000000000, 5.0000000000000000 15.0000000000000000))"
+                      "POLYGON ((5.0000000000000000 15.0000000000000000, 5.0000000000000000 20.0000000000000000, 5.0000000000000000 25.0000000000000000, 15.0000000000000000 25.0000000000000000, 15.0000000000000000 10.0000000000000000, 10.0000000000000000 5.0000000000000000, 5.0000000000000000 5.0000000000000000, 0.0000000000000000 5.0000000000000000, 0.0000000000000000 15.0000000000000000, 5.0000000000000000 15.0000000000000000))"
                   ));
 }
 

tests/unit/operation/buffer/BufferOpTest.cpp

@@ -95,7 +95,24 @@ void object::test<2>
     ensure_not(gBuffer->isEmpty());
     ensure(gBuffer->isValid());
     ensure_equals(gBuffer->getGeometryTypeId(), geos::geom::GEOS_POLYGON);
-    ensure(gBuffer->getNumPoints() > std::size_t(32));
+    ensure_equals(gBuffer->getNumPoints(), 33u);
+    auto coords = gBuffer->getCoordinates();
+    ensure_equals(coords->getSize(), 33u);
+
+    // Check four sides to check they are exactly on unit circle
+    auto coord = coords->getAt(0);
+    ensure_equals(coord.x, 1.0);
+    ensure_equals(coord.y, 0.0);
+    coord = coords->getAt(8);
+    ensure_equals(coord.x, 0.0);
+    ensure_equals(coord.y, -1.0);
+    coord = coords->getAt(16);
+    ensure_equals(coord.x, -1.0);
+    ensure_equals(coord.y, 0.0);
+    coord = coords->getAt(24);
+    ensure_equals(coord.x, 0.0);
+    ensure_equals(coord.y, 1.0);
+
 }
 
 template<>