This is an automated email from the ASF dual-hosted git repository.
mattjuntunen pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/commons-geometry.git
The following commit(s) were added to refs/heads/master by this push:
new c1c67f5 GEOMETRY-100: adding special cases for centroid computations
for very small areas
c1c67f5 is described below
commit c1c67f5b6cd782d3f01d383f0070ed59107f742b
Author: Matt Juntunen <mattjuntu...@apache.org>
AuthorDate: Sat Jul 11 16:05:49 2020 -0400
GEOMETRY-100: adding special cases for centroid computations for very small
areas
---
.../geometry/spherical/twod/ConvexArea2S.java | 167 ++++++++++++++++++--
.../geometry/spherical/twod/GreatCircle.java | 3 +-
.../geometry/spherical/twod/RegionBSPTree2S.java | 26 +++-
.../geometry/spherical/twod/ConvexArea2STest.java | 76 ++++++++-
.../spherical/twod/RegionBSPTree2STest.java | 172 +++++++++++++++++++--
5 files changed, 405 insertions(+), 39 deletions(-)
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.ArrayList;
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 @@ public final class ConvexArea2S extends
AbstractConvexHyperplaneBoundedRegion<Po
/** 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 @@ public final class ConvexArea2S extends
AbstractConvexHyperplaneBoundedRegion<Po
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 @@ public final class ConvexArea2S extends
AbstractConvexHyperplaneBoundedRegion<Po
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 @@ import org.apache.commons.numbers.angle.PlaneAngleRadians;
* 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 @@ public class RegionBSPTree2S extends
AbstractRegionBSPTree<Point2S, RegionBSPTre
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 @@ public class ConvexArea2STest {
}
@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 @@ public class ConvexArea2STest {
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.Arrays;
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 @@ public class RegionBSPTree2STest {
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 @@ public class RegionBSPTree2STest {
}
@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 @@ public class RegionBSPTree2STest {
}
@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 @@ public class RegionBSPTree2STest {
@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 @@ public class RegionBSPTree2STest {
{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 @@ public class RegionBSPTree2STest {
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 @@ public class RegionBSPTree2STest {
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 @@ public class RegionBSPTree2STest {
// 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 @@ public class RegionBSPTree2STest {
}
}
- 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 @@ public class RegionBSPTree2STest {
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 @@ public class RegionBSPTree2STest {
pts.add(rotate.apply(pts.get(i - 1)));
}
- return GreatArcPath.fromVertexLoop(pts, TEST_PRECISION).toTree();
+ return GreatArcPath.fromVertexLoop(pts, precision).toTree();
}
}
