From e4a6e1af65f7d60fa63d7666f6cb000e365bb50f Mon Sep 17 00:00:00 2001 From: Mike Taves Date: Wed, 8 Nov 2023 10:05:13 +1300 Subject: [PATCH] Add Angle::SinCos to avoid small errors, e.g. with buffer operations (#978) --- NEWS.md | 1 + include/geos/algorithm/Angle.h | 21 ++++++++++ src/algorithm/Angle.cpp | 2 +- src/operation/buffer/BufferParameters.cpp | 3 +- .../buffer/OffsetSegmentGenerator.cpp | 34 +++++++++------- src/util/GeometricShapeFactory.cpp | 31 ++++++++------- tests/unit/algorithm/AngleTest.cpp | 39 +++++++++++++++++++ tests/unit/capi/GEOSBufferTest.cpp | 8 ++-- tests/unit/operation/buffer/BufferOpTest.cpp | 19 ++++++++- 9 files changed, 123 insertions(+), 35 deletions(-) diff --git a/NEWS.md b/NEWS.md index 08752b6b11..a77bd1d326 100644 --- a/NEWS.md +++ b/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 diff --git a/include/geos/algorithm/Angle.h b/include/geos/algorithm/Angle.h index ddcdedb7b7..2c1a5edaba 100644 --- a/include/geos/algorithm/Angle.h +++ b/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; + } + }; diff --git a/src/algorithm/Angle.cpp b/src/algorithm/Angle.cpp index 3b0aebf41e..7069607f55 100644 --- a/src/algorithm/Angle.cpp +++ b/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; diff --git a/src/operation/buffer/BufferParameters.cpp b/src/operation/buffer/BufferParameters.cpp index 8c97a51c30..76c6c13e91 100644 --- a/src/operation/buffer/BufferParameters.cpp +++ b/src/operation/buffer/BufferParameters.cpp @@ -19,6 +19,7 @@ #include // for std::abs() #include // for cos +#include #include #include @@ -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); } diff --git a/src/operation/buffer/OffsetSegmentGenerator.cpp b/src/operation/buffer/OffsetSegmentGenerator.cpp index 729f6c7ac0..db7061b276 100644 --- a/src/operation/buffer/OffsetSegmentGenerator.cpp +++ b/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); } diff --git a/src/util/GeometricShapeFactory.cpp b/src/util/GeometricShapeFactory.cpp index bfe485460d..f4d714650b 100644 --- a/src/util/GeometricShapeFactory.cpp +++ b/src/util/GeometricShapeFactory.cpp @@ -17,6 +17,7 @@ * **********************************************************************/ +#include #include #include #include @@ -31,6 +32,7 @@ #include +using namespace geos::algorithm; using namespace geos::geom; namespace geos { @@ -134,10 +136,11 @@ GeometricShapeFactory::createCircle() auto pts = detail::make_unique(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(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(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); diff --git a/tests/unit/algorithm/AngleTest.cpp b/tests/unit/algorithm/AngleTest.cpp index 86a942b6c7..372b8ebeb0 100644 --- a/tests/unit/algorithm/AngleTest.cpp +++ b/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 diff --git a/tests/unit/capi/GEOSBufferTest.cpp b/tests/unit/capi/GEOSBufferTest.cpp index 0c7cab9629..f47a2237f9 100644 --- a/tests/unit/capi/GEOSBufferTest.cpp +++ b/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))" )); } diff --git a/tests/unit/operation/buffer/BufferOpTest.cpp b/tests/unit/operation/buffer/BufferOpTest.cpp index 4d4c693d94..a324bfb2ad 100644 --- a/tests/unit/operation/buffer/BufferOpTest.cpp +++ b/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<>