Skip to content

Commit

Permalink
Add Angle::SinCos to find exact angles on unit circle, avoid small er…
Browse files Browse the repository at this point in the history 
…rors
  • Loading branch information
@mwtoews
mwtoews committed Oct 30, 2023
1 parent d7a5a51 commit 3f9e832
Show file tree
Hide file tree
Showing 8 changed files with 183 additions and 35 deletions.
13 changes: 13 additions & 0 deletions include/geos/algorithm/Angle.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -216,6 +216,19 @@ 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(M_PI)
/// and std::cos(M_PI_2), this method will return exactly zero.
///
/// @param ang the input angle (in radians)
/// @param rSin the result of sin(ang)
/// @param rCos the result of cos(ang)
/// @return always 0
///
static int SinCos(const double ang, double& rSin, double& rCos);
};


Expand Down
67 changes: 66 additions & 1 deletion src/algorithm/Angle.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -195,12 +195,77 @@ Angle::diff(double ang1, double ang2)
}

if(delAngle > MATH_PI) {
delAngle = (2 * MATH_PI) - delAngle;
delAngle = PI_TIMES_2 - delAngle;
}

return delAngle;
}

/* public static */
int
Angle::SinCos(const double ang, double& rSin, double& rCos)
{
// look-up exact angles very close to 45 degree angles on unit circle
if( std::fabs(std::fmod(ang, PI_OVER_4)) < 4e-16) {
int quadrant = static_cast<int>(std::floor(ang / PI_OVER_2)) % 4;
if (quadrant < 0)
quadrant += 4;

if( std::fabs(std::fmod(ang, PI_OVER_2)) < 4e-16) {
// very close to 90 degree angles on unit circle
switch (quadrant) {
case 0:
rSin = 0.0;
rCos = 1.0;
break;
case 1:
rSin = 1.0;
rCos = 0.0;
break;
case 2:
rSin = 0.0;
rCos = -1.0;
break;
case 3:
rSin = -1.0;
rCos = 0.0;
break;
default:
std::cerr << "doom!";
}
} else {
// remainder are very close to 45 degree angles
const double SQRT2_OVER_2 = std::sqrt(2.0) / 2.0;
switch (quadrant) {
case 0:
rSin = +SQRT2_OVER_2;
rCos = +SQRT2_OVER_2;
break;
case 1:
rSin = +SQRT2_OVER_2;
rCos = -SQRT2_OVER_2;
break;
case 2:
rSin = -SQRT2_OVER_2;
rCos = -SQRT2_OVER_2;
break;
case 3:
rSin = -SQRT2_OVER_2;
rCos = +SQRT2_OVER_2;
break;
default:
std::cerr << "doom!";
}
}
} else {
// calculate other angles; may be optimized with FSINCOS instruction
rSin = std::sin(ang);
rCos = std::cos(ang);
}

return 0;
}

} // namespace geos.algorithm
} //namespace geos

3 changes: 2 additions & 1 deletion src/operation/buffer/BufferParameters.cpp
View file
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
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
43 changes: 43 additions & 0 deletions tests/unit/algorithm/AngleTest.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -162,6 +162,49 @@ void object::test<5>
TOL);
}

// testSinCos()
template<>
template<>
void object::test<6>
()
{
double rSin, rCos;
int res;

// -720 to 720 degrees with 1 degree increments
for (int angdeg = -720; angdeg <= 720; angdeg++) {
double ang = Angle::toRadians(angdeg);

res = Angle::SinCos(ang, rSin, rCos);
ensure_equals(res, 0);

double cSin, cCos;
cSin = std::sin(ang);
cCos = std::cos(ang);
// zap near-zero results to zero
if (std::fabs(cSin) < 5e-16)
cSin = 0.0;
if (std::fabs(cCos) < 5e-16)
cCos = 0.0;

ensure_equals(angdeg, rSin, cSin);
ensure_equals(angdeg, rCos, cCos);

}

// use radian increments that don't snap to special angles
for (double angrad = -6.3; angrad < 6.3; angrad += 0.013) {

res = Angle::SinCos(angrad, rSin, rCos);
ensure_equals(res, 0);

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
Original file line number Diff line number Diff line change
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
Loading

0 comments on commit 3f9e832

Please sign in to comment.