Merge remote-tracking branch 'origin/master' into 1.0-beta1-release; including changes made to master since rc2
diff --git a/commons-geometry-core/src/main/java/org/apache/commons/geometry/core/internal/GeometryInternalException.java b/commons-geometry-core/src/main/java/org/apache/commons/geometry/core/internal/GeometryInternalError.java
similarity index 93%
rename from commons-geometry-core/src/main/java/org/apache/commons/geometry/core/internal/GeometryInternalException.java
rename to commons-geometry-core/src/main/java/org/apache/commons/geometry/core/internal/GeometryInternalError.java
index caf209e..bc52f62 100644
--- a/commons-geometry-core/src/main/java/org/apache/commons/geometry/core/internal/GeometryInternalException.java
+++ b/commons-geometry-core/src/main/java/org/apache/commons/geometry/core/internal/GeometryInternalError.java
@@ -20,7 +20,7 @@
* This class is not intended to be part of the public API and should
* never be seen by the user when algorithms are functioning correctly.
*/
-public class GeometryInternalException extends IllegalStateException {
+public class GeometryInternalError extends IllegalStateException {
/** Error message used for exceptions of this type. */
private static final String ERROR_MSG = "An internal geometry error occurred. This most often indicates an " +
@@ -32,7 +32,7 @@
/** Simple constructor with a default error message.
*/
- public GeometryInternalException() {
+ public GeometryInternalError() {
super(ERROR_MSG);
}
}
diff --git a/commons-geometry-core/src/test/java/org/apache/commons/geometry/core/internal/GeometryInternalExceptionTest.java b/commons-geometry-core/src/test/java/org/apache/commons/geometry/core/internal/GeometryInternalErrorTest.java
similarity index 90%
rename from commons-geometry-core/src/test/java/org/apache/commons/geometry/core/internal/GeometryInternalExceptionTest.java
rename to commons-geometry-core/src/test/java/org/apache/commons/geometry/core/internal/GeometryInternalErrorTest.java
index 4780365..87198e1 100644
--- a/commons-geometry-core/src/test/java/org/apache/commons/geometry/core/internal/GeometryInternalExceptionTest.java
+++ b/commons-geometry-core/src/test/java/org/apache/commons/geometry/core/internal/GeometryInternalErrorTest.java
@@ -19,12 +19,12 @@
import org.junit.Assert;
import org.junit.Test;
-public class GeometryInternalExceptionTest {
+public class GeometryInternalErrorTest {
@Test
public void testMessage() {
// act
- final GeometryInternalException err = new GeometryInternalException();
+ final GeometryInternalError err = new GeometryInternalError();
// assert
final String msg = "An internal geometry error occurred. This most often indicates an " +
diff --git a/commons-geometry-enclosing/src/main/java/org/apache/commons/geometry/enclosing/WelzlEncloser.java b/commons-geometry-enclosing/src/main/java/org/apache/commons/geometry/enclosing/WelzlEncloser.java
index d94a0aa..9f9f772 100755
--- a/commons-geometry-enclosing/src/main/java/org/apache/commons/geometry/enclosing/WelzlEncloser.java
+++ b/commons-geometry-enclosing/src/main/java/org/apache/commons/geometry/enclosing/WelzlEncloser.java
@@ -20,7 +20,7 @@
import java.util.List;
import org.apache.commons.geometry.core.Point;
-import org.apache.commons.geometry.core.internal.GeometryInternalException;
+import org.apache.commons.geometry.core.internal.GeometryInternalError;
import org.apache.commons.geometry.core.precision.DoublePrecisionContext;
/** Class implementing Emo Welzl's algorithm to find the smallest enclosing ball in linear time.
@@ -98,7 +98,7 @@
ball = moveToFrontBall(extreme, extreme.size(), support);
if (precision.lt(ball.getRadius(), savedBall.getRadius())) {
// this should never happen
- throw new GeometryInternalException();
+ throw new GeometryInternalError();
}
// it was an interesting point, move it to the front
diff --git a/commons-geometry-euclidean/src/main/java/org/apache/commons/geometry/euclidean/threed/rotation/QuaternionRotation.java b/commons-geometry-euclidean/src/main/java/org/apache/commons/geometry/euclidean/threed/rotation/QuaternionRotation.java
index 65f34c1..62e950a 100644
--- a/commons-geometry-euclidean/src/main/java/org/apache/commons/geometry/euclidean/threed/rotation/QuaternionRotation.java
+++ b/commons-geometry-euclidean/src/main/java/org/apache/commons/geometry/euclidean/threed/rotation/QuaternionRotation.java
@@ -19,7 +19,7 @@
import java.util.Objects;
import java.util.function.DoubleFunction;
-import org.apache.commons.geometry.core.internal.GeometryInternalException;
+import org.apache.commons.geometry.core.internal.GeometryInternalError;
import org.apache.commons.geometry.euclidean.internal.Vectors;
import org.apache.commons.geometry.euclidean.threed.AffineTransformMatrix3D;
import org.apache.commons.geometry.euclidean.threed.Vector3D;
@@ -382,7 +382,7 @@
}
// all possibilities should have been covered above
- throw new GeometryInternalException();
+ throw new GeometryInternalError();
}
/** Get a sequence of angles around the given Tait-Bryan axes that produce a rotation equivalent
diff --git a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/ConvexArea2S.java b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/ConvexArea2S.java
index a144a66..fc9de10 100644
--- a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/ConvexArea2S.java
+++ b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/ConvexArea2S.java
@@ -20,6 +20,7 @@
import java.util.Arrays;
import java.util.Collection;
import java.util.Collections;
+import java.util.Iterator;
import java.util.List;
import java.util.stream.Collectors;
import java.util.stream.Stream;
@@ -47,6 +48,12 @@
/** Constant containing the area of half of the spherical space. */
private static final double HALF_SIZE = PlaneAngleRadians.TWO_PI;
+ /** Empirically determined threshold for computing the weighted centroid vector using the
+ * triangle fan approach. Areas with boundary sizes under this value use the triangle fan
+ * method to increase centroid accuracy.
+ */
+ private static final double TRIANGLE_FAN_CENTROID_COMPUTE_THRESHOLD = 1e-2;
+
/** Construct an instance from its boundaries. Callers are responsible for ensuring
* that the given path represents the boundary of a convex area. No validation is
* performed.
@@ -132,34 +139,35 @@
return weighted == null ? null : Point2S.from(weighted);
}
- /** Returns the weighted vector for the centroid. This vector is computed by scaling the
- * pole vector of the great circle of each boundary arc by the size of the arc and summing
- * the results. By combining the weighted centroid vectors of multiple areas, a single
- * centroid can be computed for the whole group.
+ /** Return the weighted centroid vector of the area. The returned vector points in the direction of the
+ * centroid point on the surface of the unit sphere with the length of the vector proportional to the
+ * effective mass of the area at the centroid. By adding the weighted centroid vectors of multiple
+ * convex areas, a single centroid can be computed for the combined area.
* @return weighted centroid vector.
* @see <a href="https://archive.org/details/centroidinertiat00broc">
* <em>The Centroid and Inertia Tensor for a Spherical Triangle</em> - John E. Brock</a>
*/
Vector3D getWeightedCentroidVector() {
final List<GreatArc> arcs = getBoundaries();
- switch (arcs.size()) {
+ final int numBoundaries = arcs.size();
+
+ switch (numBoundaries) {
case 0:
// full space; no centroid
return null;
case 1:
- // hemisphere; centroid is the pole of the hemisphere
- final GreatArc singleArc = arcs.get(0);
- return singleArc.getCircle().getPole().withNorm(singleArc.getSize());
+ // hemisphere
+ return computeHemisphereWeightedCentroidVector(arcs.get(0));
+ case 2:
+ // lune
+ return computeLuneWeightedCentroidVector(arcs.get(0), arcs.get(1));
default:
- // 2 or more sides; use an extension of the approach outlined here:
- // https://archive.org/details/centroidinertiat00broc
- // In short, the centroid is the sum of the pole vectors of each side
- // multiplied by their arc lengths.
- Vector3D centroid = Vector3D.ZERO;
- for (final GreatArc arc : getBoundaries()) {
- centroid = centroid.add(arc.getCircle().getPole().withNorm(arc.getSize()));
+ // triangle or other convex polygon
+ if (getBoundarySize() < TRIANGLE_FAN_CENTROID_COMPUTE_THRESHOLD) {
+ return computeTriangleFanWeightedCentroidVector(arcs);
}
- return centroid;
+
+ return computeArcPoleWeightedCentroidVector(arcs);
}
}
@@ -313,4 +321,131 @@
full() :
new ConvexArea2S(arcs);
}
+
+ /** Compute the weighted centroid vector for the hemisphere formed by the given arc.
+ * @param arc arc defining the hemisphere
+ * @return the weighted centroid vector for the hemisphere
+ * @see #getWeightedCentroidVector()
+ */
+ private static Vector3D computeHemisphereWeightedCentroidVector(final GreatArc arc) {
+ return arc.getCircle().getPole().withNorm(HALF_SIZE);
+ }
+
+ /** Compute the weighted centroid vector for the lune formed by the given arcs.
+ * @param a first arc for the lune
+ * @param b second arc for the lune
+ * @return the weighted centroid vector for the lune
+ * @see #getWeightedCentroidVector()
+ */
+ private static Vector3D computeLuneWeightedCentroidVector(final GreatArc a, final GreatArc b) {
+ final Point2S aMid = a.getCentroid();
+ final Point2S bMid = b.getCentroid();
+
+ // compute the centroid vector as the exact center of the lune to avoid inaccurate
+ // results with very small regions
+ final Vector3D.Unit centroid = aMid.slerp(bMid, 0.5).getVector();
+
+ // compute the weight using the reverse of the algorithm from computeArcPoleWeightedCentroidVector()
+ final double weight =
+ (a.getSize() * centroid.dot(a.getCircle().getPole())) +
+ (b.getSize() * centroid.dot(b.getCircle().getPole()));
+
+ return centroid.withNorm(weight);
+ }
+
+ /** Compute the weighted centroid vector for the triangle or polygon formed by the given arcs
+ * by adding together the arc pole vectors multiplied by their respective arc lengths. This
+ * algorithm is described in the paper <a href="https://archive.org/details/centroidinertiat00broc">
+ * <em>The Centroid and Inertia Tensor for a Spherical Triangle</em></a> by John E Brock.
+ *
+ * <p>Note: This algorithm works well in general but is susceptible to floating point errors
+ * on very small areas. In these cases, the computed centroid may not be in the expected location
+ * and may even be outside of the area. The {@link #computeTriangleFanWeightedCentroidVector(List)}
+ * method can produce more accurate results in these cases.</p>
+ * @param arcs boundary arcs for the area
+ * @return the weighted centroid vector for the area
+ * @see #computeTriangleFanWeightedCentroidVector(List)
+ */
+ private static Vector3D computeArcPoleWeightedCentroidVector(final List<GreatArc> arcs) {
+ Vector3D centroid = Vector3D.ZERO;
+
+ Vector3D arcContribution;
+ for (final GreatArc arc : arcs) {
+ arcContribution = arc.getCircle().getPole().withNorm(arc.getSize());
+
+ centroid = centroid.add(arcContribution);
+ }
+
+ return centroid;
+ }
+
+ /** Compute the weighted centroid vector for the triangle or polygon formed by the given arcs
+ * using a triangle fan approach. This method is specifically designed for use with areas of very small size,
+ * where use of the standard algorithm from {@link ##computeArcPoleWeightedCentroidVector(List))} can produce
+ * inaccurate results. The algorithm proceeds as follows:
+ * <ol>
+ * <li>The polygon is divided into spherical triangles using a triangle fan.</li>
+ * <li>For each triangle, the vectors of the 3 spherical points are added together to approximate the direction
+ * of the spherical centroid. This ensures that the computed centroid lies within the area.</li>
+ * <li>The length of the weighted centroid vector is determined by computing the sum of the contributions that
+ * each arc in the triangle would make to the centroid using the algorithm from
+ * {@link ##computeArcPoleWeightedCentroidVector(List)}. This essentially performs part of that algorithm in
+ * reverse: given a centroid direction, compute the contribution that each arc makes.</li>
+ * <li>The sum of the weighted centroid vectors for each triangle is computed and returned.</li>
+ * </ol>
+ * @param arcs boundary arcs for the area; must contain at least 3 arcs
+ * @return the weighted centroid vector for the area
+ * @see ##computeArcPoleWeightedCentroidVector(List)
+ */
+ private static Vector3D computeTriangleFanWeightedCentroidVector(final List<GreatArc> arcs) {
+ final Iterator<GreatArc> arcIt = arcs.iterator();
+
+ final Point2S p0 = arcIt.next().getStartPoint();
+ final Vector3D.Unit v0 = p0.getVector();
+
+ Vector3D areaCentroid = Vector3D.ZERO;
+
+ GreatArc arc;
+ Point2S p1;
+ Point2S p2;
+ Vector3D.Unit v1;
+ Vector3D.Unit v2;
+ Vector3D.Unit triangleCentroid;
+ double triangleCentroidLen;
+ while (arcIt.hasNext()) {
+ arc = arcIt.next();
+
+ if (!arc.contains(p0)) {
+ p1 = arc.getStartPoint();
+ p2 = arc.getEndPoint();
+
+ v1 = p1.getVector();
+ v2 = p2.getVector();
+
+ triangleCentroid = v0.add(v1).add(v2).normalize();
+ triangleCentroidLen =
+ computeArcCentroidContribution(v0, v1, triangleCentroid) +
+ computeArcCentroidContribution(v1, v2, triangleCentroid) +
+ computeArcCentroidContribution(v2, v0, triangleCentroid);
+
+ areaCentroid = areaCentroid.add(triangleCentroid.withNorm(triangleCentroidLen));
+ }
+ }
+
+ return areaCentroid;
+ }
+
+ /** Compute the contribution made by a single arc to a weighted centroid vector.
+ * @param a first point in the arc
+ * @param b second point in the arc
+ * @param triangleCentroid the centroid vector for the area
+ * @return the contribution made by the arc {@code ab} to the length of the weighted centroid vector
+ */
+ private static double computeArcCentroidContribution(final Vector3D.Unit a, final Vector3D.Unit b,
+ final Vector3D.Unit triangleCentroid) {
+ final double arcLength = a.angle(b);
+ final Vector3D.Unit planeNormal = a.cross(b).normalize();
+
+ return arcLength * triangleCentroid.dot(planeNormal);
+ }
}
diff --git a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/GreatCircle.java b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/GreatCircle.java
index c1e830c..4847449 100644
--- a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/GreatCircle.java
+++ b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/GreatCircle.java
@@ -32,7 +32,8 @@
* intersection of a sphere with a plane that passes through its center. It is
* the largest diameter circle that can be drawn on the sphere and partitions the
* sphere into two hemispheres. The vectors {@code u} and {@code v} lie in the great
- * circle plane, while the vector {@code w} (the pole) is perpendicular to it.
+ * circle plane, while the vector {@code w} (the pole) is perpendicular to it. The
+ * pole vector points toward the <em>minus</em> side of the hyperplane.
*
* <p>Instances of this class are guaranteed to be immutable.</p>
* @see GreatCircles
diff --git a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2S.java b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2S.java
index 5a9d20e..2b8cf2e 100644
--- a/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2S.java
+++ b/commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2S.java
@@ -171,16 +171,32 @@
final DoublePrecisionContext precision = ((GreatArc) getRoot().getCut()).getPrecision();
double sizeSum = 0;
- Vector3D centroidVector = Vector3D.ZERO;
+ Vector3D centroidVectorSum = Vector3D.ZERO;
+
+ Vector3D centroidVector;
+ double maxCentroidVectorWeightSq = 0.0;
for (final ConvexArea2S area : areas) {
sizeSum += area.getSize();
- centroidVector = centroidVector.add(area.getWeightedCentroidVector());
+
+ centroidVector = area.getWeightedCentroidVector();
+ maxCentroidVectorWeightSq = Math.max(maxCentroidVectorWeightSq, centroidVector.normSq());
+
+ centroidVectorSum = centroidVectorSum.add(centroidVector);
}
- final Point2S centroid = centroidVector.eq(Vector3D.ZERO, precision) ?
- null :
- Point2S.from(centroidVector);
+ // Convert the weighted centroid vector to a point on the sphere surface. If the centroid vector
+ // length is less than the max length of the combined convex areas and the vector itself is
+ // equivalent to zero, then we know that there are opposing and approximately equal areas in the
+ // region, resulting in an indeterminate centroid. This would occur, for example, if there were
+ // equal areas around each pole.
+ final Point2S centroid;
+ if (centroidVectorSum.normSq() < maxCentroidVectorWeightSq &&
+ centroidVectorSum.eq(Vector3D.ZERO, precision)) {
+ centroid = null;
+ } else {
+ centroid = Point2S.from(centroidVectorSum);
+ }
return new RegionSizeProperties<>(sizeSum, centroid);
}
diff --git a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/ConvexArea2STest.java b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/ConvexArea2STest.java
index 6683be5..ddcb505 100644
--- a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/ConvexArea2STest.java
+++ b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/ConvexArea2STest.java
@@ -538,6 +538,76 @@
}
@Test
+ public void testGetCentroid_diminishingLunes() {
+ // arrange
+ final double eps = 1e-14;
+ final DoublePrecisionContext precision = new EpsilonDoublePrecisionContext(eps);
+
+ final double centerAz = 1;
+ final double centerPolar = 0.5 * Math.PI;
+ final Point2S center = Point2S.of(centerAz, centerPolar);
+ final Point2S pole = Point2S.PLUS_K;
+
+ final double startOffset = PlaneAngleRadians.PI_OVER_TWO;
+ final double minOffset = 1e-14;
+
+ ConvexArea2S area;
+ Point2S p1;
+ Point2S p2;
+ Point2S centroid;
+ for (double offset = startOffset; offset > minOffset; offset *= 0.5) {
+ p1 = Point2S.of(centerAz - offset, centerPolar);
+ p2 = Point2S.of(centerAz + offset, centerPolar);
+
+ area = ConvexArea2S.fromBounds(
+ GreatCircles.fromPoints(pole, p1, precision),
+ GreatCircles.fromPoints(p2, pole, precision));
+
+ // act
+ centroid = area.getCentroid();
+
+ // assert
+ Assert.assertTrue(area.contains(centroid));
+ SphericalTestUtils.assertPointsEq(center, centroid, TEST_EPS);
+ }
+ }
+
+ @Test
+ public void testGetCentroid_diminishingSquares() {
+ // arrange
+ final double eps = 1e-14;
+ final DoublePrecisionContext precision = new EpsilonDoublePrecisionContext(eps);
+
+ final double centerAz = 1;
+ final double centerPolar = 0.5 * Math.PI;
+ final Point2S center = Point2S.of(centerAz, centerPolar);
+
+ final double minOffset = 1e-14;
+
+ ConvexArea2S area;
+ Point2S p1;
+ Point2S p2;
+ Point2S p3;
+ Point2S p4;
+ Point2S centroid;
+ for (double offset = 0.5; offset > minOffset; offset *= 0.5) {
+ p1 = Point2S.of(centerAz, centerPolar - offset);
+ p2 = Point2S.of(centerAz - offset, centerPolar);
+ p3 = Point2S.of(centerAz, centerPolar + offset);
+ p4 = Point2S.of(centerAz + offset, centerPolar);
+
+ area = ConvexArea2S.fromVertexLoop(Arrays.asList(p1, p2, p3, p4), precision);
+
+ // act
+ centroid = area.getCentroid();
+
+ // assert
+ Assert.assertTrue(area.contains(centroid));
+ SphericalTestUtils.assertPointsEq(center, centroid, TEST_EPS);
+ }
+ }
+
+ @Test
public void testBoundaryStream() {
// arrange
final GreatCircle circle = GreatCircles.fromPole(Vector3D.Unit.PLUS_X, TEST_PRECISION);
@@ -823,8 +893,10 @@
final ConvexArea2S plus = split.getPlus();
final double plusSize = plus.getSize();
- final Point2S computedCentroid = Point2S.from(minus.getWeightedCentroidVector()
- .add(plus.getWeightedCentroidVector()));
+ final Vector3D minusWeightedCentroid = minus.getWeightedCentroidVector();
+ final Vector3D plusWeightedCentroid = plus.getWeightedCentroidVector();
+
+ final Point2S computedCentroid = Point2S.from(minusWeightedCentroid.add(plusWeightedCentroid));
Assert.assertEquals(size, minusSize + plusSize, TEST_EPS);
SphericalTestUtils.assertPointsEq(centroid, computedCentroid, TEST_EPS);
diff --git a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2STest.java b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2STest.java
index fb0784c..1b8032f 100644
--- a/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2STest.java
+++ b/commons-geometry-spherical/src/test/java/org/apache/commons/geometry/spherical/twod/RegionBSPTree2STest.java
@@ -21,6 +21,7 @@
import java.util.Collections;
import java.util.List;
import java.util.stream.Collectors;
+import java.util.stream.Stream;
import org.apache.commons.geometry.core.GeometryTestUtils;
import org.apache.commons.geometry.core.RegionLocation;
@@ -361,7 +362,11 @@
SphericalTestUtils.assertPointsEq(Point2S.MINUS_I,
tree.project(Point2S.of(PlaneAngleRadians.PI + 0.5, PlaneAngleRadians.PI_OVER_TWO)), TEST_EPS);
- SphericalTestUtils.assertPointsEq(Point2S.PLUS_K, tree.project(tree.getCentroid()), TEST_EPS);
+ final Point2S centroid = tree.getCentroid();
+ SphericalTestUtils.assertPointsEq(Point2S.PLUS_K,
+ tree.project(centroid.slerp(Point2S.PLUS_K, 1e-10)), TEST_EPS);
+ SphericalTestUtils.assertPointsEq(Point2S.PLUS_J,
+ tree.project(centroid.slerp(Point2S.PLUS_J, 1e-10)), TEST_EPS);
}
@Test
@@ -560,6 +565,42 @@
}
@Test
+ public void testGeometricProperties_polygonWithHole_small() {
+ // arrange
+ final Point2S center = Point2S.of(0.5, 2);
+
+ final double outerRadius = 1e-5;
+ final double innerRadius = 1e-7;
+
+ final RegionBSPTree2S outer = buildDiamond(center, outerRadius);
+ final RegionBSPTree2S inner = buildDiamond(center, innerRadius);
+
+ // rotate the inner diamond a quarter turn to become a square
+ inner.transform(Transform2S.createRotation(center, 0.25 * Math.PI));
+
+ // act
+ final RegionBSPTree2S tree = RegionBSPTree2S.empty();
+ tree.difference(outer, inner);
+
+ // assert
+
+ // use Euclidean approximations of the area and boundary size since those will be more accurate
+ // at these sizes
+ final double area = (2 * outerRadius * outerRadius) - (2 * innerRadius * innerRadius);
+ Assert.assertEquals(area, tree.getSize(), TEST_EPS);
+
+ final double outerSideLength = Math.hypot(outerRadius, outerRadius);
+ final double innerSideLength = Math.hypot(innerRadius, innerRadius);
+ final double boundarySize = 4 * (outerSideLength + innerSideLength);
+ Assert.assertEquals(boundarySize, tree.getBoundarySize(), TEST_EPS);
+
+ SphericalTestUtils.assertPointsEq(center, tree.getCentroid(), TEST_EPS);
+ checkCentroidConsistency(tree);
+
+ SphericalTestUtils.checkClassify(tree, RegionLocation.OUTSIDE, center);
+ }
+
+ @Test
public void testGeometricProperties_polygonWithHole_complex() {
// arrange
final Point2S center = Point2S.of(0.5, 2);
@@ -602,6 +643,59 @@
}
@Test
+ public void testGeometricProperties_smallRightTriangle() {
+ // arrange
+ final double azOffset = 1e-5;
+ final double polarOffset = 1e-6;
+
+ final double minAz = 0;
+ final double maxAz = minAz + azOffset;
+ final double maxPolar = PlaneAngleRadians.PI_OVER_TWO;
+ final double minPolar = maxPolar - polarOffset;
+
+ final Point2S p0 = Point2S.of(minAz, maxPolar);
+ final Point2S p1 = Point2S.of(maxAz, maxPolar);
+ final Point2S p2 = Point2S.of(maxAz, minPolar);
+
+ // act
+ final RegionBSPTree2S tree = GreatArcPath.fromVertexLoop(Arrays.asList(p0, p1, p2), TEST_PRECISION)
+ .toTree();
+
+ // assert
+
+ // use Euclidean approximations of the area and boundary size since those will be more accurate
+ // at these sizes
+ final double expectedArea = 0.5 * azOffset * polarOffset;
+ Assert.assertEquals(expectedArea, tree.getSize(), TEST_EPS);
+
+ final double expectedBoundarySize = azOffset + polarOffset + Math.hypot(azOffset, polarOffset);
+ Assert.assertEquals(expectedBoundarySize, tree.getBoundarySize(), TEST_EPS);
+
+ Assert.assertTrue(tree.contains(tree.getCentroid()));
+ checkCentroidConsistency(tree);
+
+ SphericalTestUtils.checkClassify(tree, RegionLocation.INSIDE,
+ tree.getCentroid(),
+ Point2S.of(minAz + (0.75 * azOffset), minPolar + (0.75 * polarOffset)));
+
+ SphericalTestUtils.checkClassify(tree, RegionLocation.BOUNDARY,
+ p0, p1, p2, p0.slerp(p1, 0.5), p1.slerp(p2, 0.5), p2.slerp(p0, 0.5));
+
+ final double midAz = minAz + (0.5 * azOffset);
+ final double pastMinAz = minAz - azOffset;
+ final double pastMaxAz = maxAz + azOffset;
+
+ final double midPolar = minPolar + (0.5 * polarOffset);
+ final double pastMinPolar = minPolar - polarOffset;
+ final double pastMaxPolar = maxPolar + polarOffset;
+
+ SphericalTestUtils.checkClassify(tree, RegionLocation.OUTSIDE,
+ tree.getCentroid().antipodal(),
+ Point2S.of(pastMinAz, midPolar), Point2S.of(pastMaxAz, midPolar),
+ Point2S.of(midAz, pastMinPolar), Point2S.of(midAz, pastMaxPolar));
+ }
+
+ @Test
public void testGeometricProperties_equalAndOppositeRegions() {
// arrange
final Point2S center = Point2S.PLUS_I;
@@ -761,7 +855,7 @@
@Test
public void testGeographicMap() {
// arrange
- final RegionBSPTree2S continental = latLongToTree(new double[][] {
+ final RegionBSPTree2S continental = latLongToTree(TEST_PRECISION, new double[][] {
{51.14850, 2.51357}, {50.94660, 1.63900}, {50.12717, 1.33876}, {49.34737, -0.98946},
{49.77634, -1.93349}, {48.64442, -1.61651}, {48.90169, -3.29581}, {48.68416, -4.59234},
{47.95495, -4.49155}, {47.57032, -2.96327}, {46.01491, -1.19379}, {44.02261, -1.38422},
@@ -771,7 +865,7 @@
{46.27298, 6.02260}, {46.72577, 6.03738}, {47.62058, 7.46675}, {49.01778, 8.09927},
{49.20195, 6.65822}, {49.44266, 5.89775}, {49.98537, 4.79922}
});
- final RegionBSPTree2S corsica = latLongToTree(new double[][] {
+ final RegionBSPTree2S corsica = latLongToTree(TEST_PRECISION, new double[][] {
{42.15249, 9.56001}, {43.00998, 9.39000}, {42.62812, 8.74600}, {42.25651, 8.54421},
{41.58361, 8.77572}, {41.38000, 9.22975}
});
@@ -797,11 +891,12 @@
final int numPts = 200;
// counterclockwise
- final RegionBSPTree2S ccw = circleToPolygon(center, radius, numPts, false);
- SphericalTestUtils.assertPointsEq(center, ccw.getCentroid(), CENTROID_EPS);
+ final RegionBSPTree2S ccw = circleToPolygon(center, radius, numPts, false, TEST_PRECISION);
+ SphericalTestUtils.assertPointsEq(center, ccw.getCentroid(), TEST_EPS);
// clockwise; centroid should just be antipodal for the circle center
- final RegionBSPTree2S cw = circleToPolygon(center, radius, numPts, true);
+ final RegionBSPTree2S cw = circleToPolygon(center, radius, numPts, true, TEST_PRECISION);
+
SphericalTestUtils.assertPointsEq(center.antipodal(), cw.getCentroid(), CENTROID_EPS);
}
@@ -815,10 +910,10 @@
final double ccwArea = 4.0 * PlaneAngleRadians.PI * Math.pow(Math.sin(radius / 2.0), 2.0);
final double cwArea = 4.0 * PlaneAngleRadians.PI - ccwArea;
- final RegionBSPTree2S ccw = circleToPolygon(center, radius, numPts, false);
+ final RegionBSPTree2S ccw = circleToPolygon(center, radius, numPts, false, TEST_PRECISION);
Assert.assertEquals("Counterclockwise size", ccwArea, ccw.getSize(), TEST_EPS);
- final RegionBSPTree2S cw = circleToPolygon(center, radius, numPts, true);
+ final RegionBSPTree2S cw = circleToPolygon(center, radius, numPts, true, TEST_PRECISION);
Assert.assertEquals("Clockwise size", cwArea, cw.getSize(), TEST_EPS);
}
@@ -831,13 +926,59 @@
// boundary size is independent from winding
final double boundary = PlaneAngleRadians.TWO_PI * Math.sin(radius);
- final RegionBSPTree2S ccw = circleToPolygon(center, radius, numPts, false);
+ final RegionBSPTree2S ccw = circleToPolygon(center, radius, numPts, false, TEST_PRECISION);
Assert.assertEquals("Counterclockwise boundary size", boundary, ccw.getBoundarySize(), 1.0e-7);
- final RegionBSPTree2S cw = circleToPolygon(center, radius, numPts, true);
+ final RegionBSPTree2S cw = circleToPolygon(center, radius, numPts, true, TEST_PRECISION);
Assert.assertEquals("Clockwise boundary size", boundary, cw.getBoundarySize(), 1.0e-7);
}
+ @Test
+ public void testSmallCircleToPolygon() {
+ // arrange
+ final double radius = 5.0e-8;
+ final Point2S center = Point2S.of(0.5, 1.5);
+ final int numPts = 100;
+
+ // act
+ final RegionBSPTree2S circle = circleToPolygon(center, radius, numPts, false, TEST_PRECISION);
+
+ // assert
+ // https://en.wikipedia.org/wiki/Spherical_cap
+ final double area = 4.0 * PlaneAngleRadians.PI * Math.pow(Math.sin(radius / 2.0), 2.0);
+ final double boundary = PlaneAngleRadians.TWO_PI * Math.sin(radius);
+
+ SphericalTestUtils.assertPointsEq(center, circle.getCentroid(), TEST_EPS);
+ Assert.assertEquals(area, circle.getSize(), TEST_EPS);
+ Assert.assertEquals(boundary, circle.getBoundarySize(), TEST_EPS);
+ }
+
+ @Test
+ public void testSmallGeographicalRectangle() {
+ // arrange
+ final double[][] vertices = {
+ {42.656216727628696, -70.61919768884546},
+ {42.65612858998112, -70.61938607250165},
+ {42.65579098923594, -70.61909615581666},
+ {42.655879126692355, -70.61890777301083}
+ };
+
+ // act
+ final RegionBSPTree2S rectangle = latLongToTree(TEST_PRECISION, vertices);
+
+ // assert
+ // approximate the centroid as average of vertices
+ final double avgLat = Stream.of(vertices).mapToDouble(v -> v[0]).average().getAsDouble();
+ final double avgLon = Stream.of(vertices).mapToDouble(v -> v[1]).average().getAsDouble();
+ final Point2S expectedCentroid = latLongToPoint(avgLat, avgLon);
+
+ SphericalTestUtils.assertPointsEq(expectedCentroid, rectangle.getCentroid(), TEST_EPS);
+
+ // expected results computed using GeographicLib (https://geographiclib.sourceforge.io/)
+ Assert.assertEquals(1.997213869978027E-11, rectangle.getSize(), TEST_EPS);
+ Assert.assertEquals(1.9669710464585642E-5, rectangle.getBoundarySize(), TEST_EPS);
+ }
+
/**
* Insert hyperplane convex subsets defining the positive quadrant area.
* @param tree
@@ -875,8 +1016,8 @@
}
}
- private static RegionBSPTree2S latLongToTree(final double[][] points) {
- final GreatArcPath.Builder pathBuilder = GreatArcPath.builder(TEST_PRECISION);
+ private static RegionBSPTree2S latLongToTree(final DoublePrecisionContext precision, final double[][] points) {
+ final GreatArcPath.Builder pathBuilder = GreatArcPath.builder(precision);
for (int i = 0; i < points.length; ++i) {
pathBuilder.append(latLongToPoint(points[i][0], points[i][1]));
@@ -972,10 +1113,11 @@
final double angleB = Math.asin(Math.sin(b) / sinC);
// use Girard's theorem
- return angleA + angleB - (0.5 * PlaneAngleRadians.PI);
+ return angleA + angleB - PlaneAngleRadians.PI_OVER_TWO;
}
- private static RegionBSPTree2S circleToPolygon(final Point2S center, final double radius, final int numPts, final boolean clockwise) {
+ private static RegionBSPTree2S circleToPolygon(final Point2S center, final double radius, final int numPts,
+ final boolean clockwise, final DoublePrecisionContext precision) {
final List<Point2S> pts = new ArrayList<>(numPts);
// get an arbitrary point on the circle boundary
@@ -990,6 +1132,6 @@
pts.add(rotate.apply(pts.get(i - 1)));
}
- return GreatArcPath.fromVertexLoop(pts, TEST_PRECISION).toTree();
+ return GreatArcPath.fromVertexLoop(pts, precision).toTree();
}
}
