From 5f83405c118b54ece0caee6b9db286c55dd0dbcc Mon Sep 17 00:00:00 2001 From: Mike Taves Date: Wed, 8 Nov 2023 12:31:19 +1300 Subject: [PATCH] Add Angle.sinSnap and .cosSnap to avoid small errors, e.g. with buffer operations --- .../org/locationtech/jts/algorithm/Angle.java | 32 ++++++++++++++--- .../operation/buffer/BufferParameters.java | 4 ++- .../buffer/OffsetSegmentGenerator.java | 24 ++++++------- .../jts/util/GeometricShapeFactory.java | 21 +++++------ .../locationtech/jts/algorithm/AngleTest.java | 36 ++++++++++++++++++- .../jts/operation/buffer/BufferTest.java | 19 ++++++++++ 6 files changed, 108 insertions(+), 28 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/Angle.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/Angle.java index 1da839b395..58dd05b34e 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/Angle.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/Angle.java @@ -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. * @@ -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); } } diff --git a/modules/core/src/main/java/org/locationtech/jts/operation/buffer/BufferParameters.java b/modules/core/src/main/java/org/locationtech/jts/operation/buffer/BufferParameters.java index e9cc6f8799..f4799fc26b 100644 --- a/modules/core/src/main/java/org/locationtech/jts/operation/buffer/BufferParameters.java +++ b/modules/core/src/main/java/org/locationtech/jts/operation/buffer/BufferParameters.java @@ -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. @@ -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); } diff --git a/modules/core/src/main/java/org/locationtech/jts/operation/buffer/OffsetSegmentGenerator.java b/modules/core/src/main/java/org/locationtech/jts/operation/buffer/OffsetSegmentGenerator.java index f1591e4451..611d5e0c9b 100644 --- a/modules/core/src/main/java/org/locationtech/jts/operation/buffer/OffsetSegmentGenerator.java +++ b/modules/core/src/main/java/org/locationtech/jts/operation/buffer/OffsetSegmentGenerator.java @@ -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 @@ -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: @@ -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, @@ -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); @@ -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); } @@ -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); @@ -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); } } @@ -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(); } diff --git a/modules/core/src/main/java/org/locationtech/jts/util/GeometricShapeFactory.java b/modules/core/src/main/java/org/locationtech/jts/util/GeometricShapeFactory.java index 7ba7dd2ec9..90f5bb26a8 100644 --- a/modules/core/src/main/java/org/locationtech/jts/util/GeometricShapeFactory.java +++ b/modules/core/src/main/java/org/locationtech/jts/util/GeometricShapeFactory.java @@ -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; @@ -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]); @@ -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); @@ -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; @@ -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); diff --git a/modules/core/src/test/java/org/locationtech/jts/algorithm/AngleTest.java b/modules/core/src/test/java/org/locationtech/jts/algorithm/AngleTest.java index 417bb1f538..9f005cef06 100644 --- a/modules/core/src/test/java/org/locationtech/jts/algorithm/AngleTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/algorithm/AngleTest.java @@ -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); } diff --git a/modules/core/src/test/java/org/locationtech/jts/operation/buffer/BufferTest.java b/modules/core/src/test/java/org/locationtech/jts/operation/buffer/BufferTest.java index 73b5509362..6e32380d8d 100644 --- a/modules/core/src/test/java/org/locationtech/jts/operation/buffer/BufferTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/operation/buffer/BufferTest.java @@ -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; @@ -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); + } + }