Skip to content

Commit

Permalink
Add Angle.sinSnap and .cosSnap to avoid small errors, e.g. with buffe…
Browse files Browse the repository at this point in the history 
…r operations (#1016)
  • Loading branch information
@mwtoews
mwtoews authored Nov 15, 2023
1 parent 84bcead commit 58d7bd5
Show file tree
Hide file tree
Showing 6 changed files with 108 additions and 28 deletions.
32 changes: 28 additions & 4 deletions modules/core/src/main/java/org/locationtech/jts/algorithm/Angle.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -313,12 +313,36 @@ public static double diff(double ang1, double ang2) {
}

if (delAngle > Math.PI) {
delAngle = (2 * Math.PI) - delAngle;
delAngle = PI_TIMES_2 - delAngle;
}

return delAngle;
}


/**
* Computes sin of an angle, snapping near-zero values to zero.
*
* @param ang the input angle (in radians)
* @return the result of the trigonometric function
*/
public static double sinSnap(double ang) {
double res = Math.sin(ang);
if (Math.abs(res) < 5e-16) return 0.0;
return res;
}

/**
* Computes cos of an angle, snapping near-zero values to zero.
*
* @param ang the input angle (in radians)
* @return the result of the trigonometric function
*/
public static double cosSnap(double ang) {
double res = Math.cos(ang);
if (Math.abs(res) < 5e-16) return 0.0;
return res;
}

/**
* Projects a point by a given angle and distance.
*
Expand All @@ -328,8 +352,8 @@ public static double diff(double ang1, double ang2) {
* @return the projected point
*/
public static Coordinate project(Coordinate p, double angle, double dist) {
double x = p.getX() + dist * Math.cos(angle);
double y = p.getY() + dist * Math.sin(angle);
double x = p.getX() + dist * Angle.cosSnap(angle);
double y = p.getY() + dist * Angle.sinSnap(angle);
return new Coordinate(x, y);
}
}
4 changes: 3 additions & 1 deletion modules/core/src/main/java/org/locationtech/jts/operation/buffer/BufferParameters.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,8 @@
*/
package org.locationtech.jts.operation.buffer;

import org.locationtech.jts.algorithm.Angle;

/**
* A value class containing the parameters which
* specify how a buffer should be constructed.
Expand Down Expand Up @@ -176,7 +178,7 @@ public void setQuadrantSegments(int quadSegs)
*/
public static double bufferDistanceError(int quadSegs)
{
double alpha = Math.PI / 2.0 / quadSegs;
double alpha = Angle.PI_OVER_2 / quadSegs;
return 1 - Math.cos(alpha / 2.0);
}

Expand Down
24 changes: 12 additions & 12 deletions modules/core/src/main/java/org/locationtech/jts/operation/buffer/OffsetSegmentGenerator.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,7 +114,7 @@ public OffsetSegmentGenerator(PrecisionModel precisionModel,

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, so don't use
Expand Down Expand Up @@ -428,7 +428,7 @@ public void addLineEndCap(Coordinate p0, 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, angle - Math.PI / 2, Orientation.CLOCKWISE, distance);
addDirectedFillet(p1, angle + Angle.PI_OVER_2, angle - Angle.PI_OVER_2, Orientation.CLOCKWISE, distance);
segList.addPt(offsetR.p1);
break;
case BufferParameters.CAP_FLAT:
Expand All @@ -439,8 +439,8 @@ public void addLineEndCap(Coordinate p0, Coordinate p1)
case BufferParameters.CAP_SQUARE:
// add a square defined by extensions of the offset segment endpoints
Coordinate squareCapSideOffset = new Coordinate();
squareCapSideOffset.x = Math.abs(distance) * Math.cos(angle);
squareCapSideOffset.y = Math.abs(distance) * Math.sin(angle);
squareCapSideOffset.x = Math.abs(distance) * Angle.cosSnap(angle);
squareCapSideOffset.y = Math.abs(distance) * Angle.sinSnap(angle);

Coordinate squareCapLOffset = new Coordinate(
offsetL.p1.x + squareCapSideOffset.x,
Expand Down Expand Up @@ -534,7 +534,7 @@ private void addLimitedMitreJoin(
Coordinate bevelMidPt = project(cornerPt, -mitreLimitDistance, dirBisector);

// direction of bevel segment (at right angle to corner bisector)
double dirBevel = Angle.normalize(dirBisector + Math.PI/2.0);
double dirBevel = Angle.normalize(dirBisector + Angle.PI_OVER_2);

// compute the candidate bevel segment by projecting both sides of the midpoint
Coordinate bevel0 = project(bevelMidPt, distance, dirBevel);
Expand Down Expand Up @@ -567,8 +567,8 @@ private void addLimitedMitreJoin(
* @return the projected point
*/
private static Coordinate project(Coordinate pt, double d, double dir) {
double x = pt.x + d * Math.cos(dir);
double y = pt.y + d * Math.sin(dir);
double x = pt.x + d * Angle.cosSnap(dir);
double y = pt.y + d * Angle.sinSnap(dir);
return new Coordinate(x, y);
}

Expand Down Expand Up @@ -607,10 +607,10 @@ private void addCornerFillet(Coordinate p, Coordinate p0, Coordinate p1, int dir
double endAngle = Math.atan2(dy1, dx1);

if (direction == Orientation.CLOCKWISE) {
if (startAngle <= endAngle) startAngle += 2.0 * Math.PI;
if (startAngle <= endAngle) startAngle += Angle.PI_TIMES_2;
}
else { // direction == COUNTERCLOCKWISE
if (startAngle >= endAngle) startAngle -= 2.0 * Math.PI;
if (startAngle >= endAngle) startAngle -= Angle.PI_TIMES_2;
}
segList.addPt(p0);
addDirectedFillet(p, startAngle, endAngle, direction, radius);
Expand Down Expand Up @@ -641,8 +641,8 @@ private void addDirectedFillet(Coordinate p, double startAngle, double endAngle,
Coordinate pt = new Coordinate();
for (int i = 0; i < nSegs; i++) {
double angle = startAngle + directionFactor * i * angleInc;
pt.x = p.x + radius * Math.cos(angle);
pt.y = p.y + radius * Math.sin(angle);
pt.x = p.x + radius * Angle.cosSnap(angle);
pt.y = p.y + radius * Angle.sinSnap(angle);
segList.addPt(pt);
}
}
Expand All @@ -655,7 +655,7 @@ public void createCircle(Coordinate p)
// add start point
Coordinate pt = new Coordinate(p.x + distance, p.y);
segList.addPt(pt);
addDirectedFillet(p, 0.0, 2.0 * Math.PI, -1, distance);
addDirectedFillet(p, 0.0, Angle.PI_TIMES_2, -1, distance);
segList.closeRing();
}

Expand Down
21 changes: 11 additions & 10 deletions modules/core/src/main/java/org/locationtech/jts/util/GeometricShapeFactory.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
*/
package org.locationtech.jts.util;

import org.locationtech.jts.algorithm.Angle;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Envelope;
import org.locationtech.jts.geom.Geometry;
Expand Down Expand Up @@ -222,8 +223,8 @@ public Polygon createEllipse()
int iPt = 0;
for (int i = 0; i < nPts; i++) {
double ang = i * (2 * Math.PI / nPts);
double x = xRadius * Math.cos(ang) + centreX;
double y = yRadius * Math.sin(ang) + centreY;
double x = xRadius * Angle.cosSnap(ang) + centreX;
double y = yRadius * Angle.sinSnap(ang) + centreY;
pts[iPt++] = coord(x, y);
}
pts[iPt] = new Coordinate(pts[0]);
Expand Down Expand Up @@ -319,16 +320,16 @@ public LineString createArc(
double centreY = env.getMinY() + yRadius;

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);

Coordinate[] pts = new Coordinate[nPts];
int iPt = 0;
for (int i = 0; i < nPts; i++) {
double ang = startAng + i * angInc;
double x = xRadius * Math.cos(ang) + centreX;
double y = yRadius * Math.sin(ang) + centreY;
double x = xRadius * Angle.cosSnap(ang) + centreX;
double y = yRadius * Angle.sinSnap(ang) + centreY;
pts[iPt++] = coord(x, y);
}
LineString line = geomFact.createLineString(pts);
Expand All @@ -353,8 +354,8 @@ public Polygon createArcPolygon(double startAng, double angExtent) {
double centreY = env.getMinY() + yRadius;

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);
// double check = angInc * nPts;
// double checkEndAng = startAng + check;
Expand All @@ -366,8 +367,8 @@ public Polygon createArcPolygon(double startAng, double angExtent) {
for (int i = 0; i < nPts; i++) {
double ang = startAng + angInc * i;

double x = xRadius * Math.cos(ang) + centreX;
double y = yRadius * Math.sin(ang) + centreY;
double x = xRadius * Angle.cosSnap(ang) + centreX;
double y = yRadius * Angle.sinSnap(ang) + centreY;
pts[iPt++] = coord(x, y);
}
pts[iPt++] = coord(centreX, centreY);
Expand Down
36 changes: 35 additions & 1 deletion modules/core/src/test/java/org/locationtech/jts/algorithm/AngleTest.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,7 +158,41 @@ public void testAngleBisector() {

assertEquals(45, Math.toDegrees(Angle.bisector(p(13,10), p(10,10), p(10,20))), 0.01);
}


public void testSinCosSnap() {

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

double rSin = Angle.sinSnap(ang);
double rCos = Angle.cosSnap(ang);

double cSin = Math.sin(ang);
double cCos = Math.cos(ang);
if ( (angdeg % 90) == 0 ) {
// not always the same for multiples of 90 degrees
assertTrue(Math.abs(rSin - cSin) < 1e-15);
assertTrue(Math.abs(rCos - cCos) < 1e-15);
} else {
assertEquals(rSin, cSin);
assertEquals(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) {

double rSin = Angle.sinSnap(angrad);
double rCos = Angle.cosSnap(angrad);

assertEquals(rSin, Math.sin(angrad));
assertEquals(rCos, Math.cos(angrad));

}
}

private static Coordinate p(double x, double y) {
return new Coordinate(x, y);
}
Expand Down
19 changes: 19 additions & 0 deletions modules/core/src/test/java/org/locationtech/jts/operation/buffer/BufferTest.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
*/
package org.locationtech.jts.operation.buffer;

import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryCollection;
import org.locationtech.jts.geom.GeometryFactory;
Expand Down Expand Up @@ -602,4 +603,22 @@ private boolean hasHole(Geometry geom) {
return false;
}

/**
* See GEOS PR https://github.com/libgeos/geos/pull/978
*/
public void testDefaultBuffer() {
Geometry g = read("POINT (0 0)").buffer(1.0);
Geometry b = g.getBoundary();
Coordinate[] coords = b.getCoordinates();
assertEquals(33, coords.length);
assertEquals(coords[0].x, 1.0);
assertEquals(coords[0].y, 0.0);
assertEquals(coords[8].x, 0.0);
assertEquals(coords[8].y, -1.0);
assertEquals(coords[16].x, -1.0);
assertEquals(coords[16].y, 0.0);
assertEquals(coords[24].x, 0.0);
assertEquals(coords[24].y, 1.0);
}

}

0 comments on commit 58d7bd5

Please sign in to comment.