Skip to content

Commit

Permalink
Add Angle::SinCos to avoid small errors, e.g. with buffer operations (#…
Browse files Browse the repository at this point in the history 
…978)
  • Loading branch information
@mwtoews
mwtoews committed Nov 7, 2023
1 parent 02f8fb1 commit e4a6e1a
Show file tree
Hide file tree
Showing 9 changed files with 123 additions and 35 deletions.
1 change: 1 addition & 0 deletions NEWS.md
View file
Expand Up @@ -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

Expand Down
21 changes: 21 additions & 0 deletions include/geos/algorithm/Angle.h
View file
Expand Up @@ -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;
}

};


Expand Down
2 changes: 1 addition & 1 deletion src/algorithm/Angle.cpp
View file
Expand Up @@ -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;
Expand Down
3 changes: 2 additions & 1 deletion src/operation/buffer/BufferParameters.cpp
View file
Expand Up @@ -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>

Expand Down Expand Up @@ -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);
}

Expand Down
34 changes: 19 additions & 15 deletions src/operation/buffer/OffsetSegmentGenerator.cpp
View file
Expand Up @@ -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,
Expand Down Expand Up @@ -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;
Expand All @@ -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,
Expand Down Expand Up @@ -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;
}
}

Expand All @@ -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);
}
}
Expand All @@ -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();
}

Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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);
}

Expand Down
31 changes: 18 additions & 13 deletions src/util/GeometricShapeFactory.cpp
View file
Expand Up @@ -17,6 +17,7 @@
*
**********************************************************************/

#include <geos/algorithm/Angle.h>
#include <geos/constants.h>
#include <geos/util/GeometricShapeFactory.h>
#include <geos/geom/Coordinate.h>
Expand All @@ -31,6 +32,7 @@
#include <memory>


using namespace geos::algorithm;
using namespace geos::geom;

namespace geos {
Expand Down Expand Up @@ -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];
Expand All @@ -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));
Expand All @@ -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);
Expand Down
39 changes: 39 additions & 0 deletions tests/unit/algorithm/AngleTest.cpp
View file
Expand Up @@ -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

8 changes: 4 additions & 4 deletions tests/unit/capi/GEOSBufferTest.cpp
View file
Expand Up @@ -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))"
));
}

Expand Down Expand Up @@ -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))"
));
}

Expand Down Expand Up @@ -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))"
));
}

Expand Down Expand Up @@ -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))"
));
}

Expand Down
19 changes: 18 additions & 1 deletion tests/unit/operation/buffer/BufferOpTest.cpp
View file
Expand Up @@ -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<>
Expand Down

0 comments on commit e4a6e1a

Please sign in to comment.